Merge remote-tracking branch 'origin/sponce_freiss_dd4hep' into rm-align-online-2

Alignment/AlignKernel/include/AlignKernel/DetectorElement.h

@@ -164,6 +164,8 @@ namespace LHCb::Alignment {
     bool operator==( const DetectorElement& other ) const { return m_data == other.m_data; }
     bool operator<( const DetectorElement& other ) const { return m_data < other.m_data; }
 
+    const auto data() const { return m_data; }
+
   private:
 #ifdef USE_DD4HEP
     ::DetectorElement m_data;

Alignment/Humboldt/options/AlignRun3Collisions.py

@@ -29,6 +29,7 @@ configure_input(options)
 
 # only configure data flow after this line !
 from Humboldt.options import usePrKalman
+#usePrKalman=False
 from Humboldt.alignment_tracking import make_align_input_besttracks
 aligntracks, alignpvs = make_align_input_besttracks(usePrKalman)
 
@@ -42,10 +43,8 @@ alignscenario = configureRun3Alignment()
 
 from Humboldt.utils import runAlignment
 from Humboldt.utils import createAlignUpdateTool, createAlignAlgorithm, getXMLWriterList
-with createAlignUpdateTool.bind(
-        logFile="alignlog_ft_modules_d0.txt"), createAlignAlgorithm.bind(
-            xmlWriters=getXMLWriterList(
-                alignscenario.SubDetectors, prefix='humb-ft-modules-d0/')):
+with createAlignAlgorithm.bind(
+        xmlWriters=getXMLWriterList(alignscenario.SubDetectors)):
     runAlignment(
         options,
         surveyConstraints=alignscenario.SurveyConstraints,

Alignment/Humboldt/python/Humboldt/alignment_tracking.py

@@ -129,10 +129,11 @@ def make_align_input_besttracks(usePrKalman=True):
             # TODO: in principle the PVs should be an optional input for the alignment
             reco = reconstruction()
     else:
+        from RecoConf.hlt2_global_reco import reconstruction as reconstruction_from_reco
+        from RecoConf.hlt2_global_reco import make_legacy_reconstruction, make_hlt2_tracks
         with reconstruction.bind(from_file=False),\
           reconstruction_from_reco.bind(make_reconstruction = make_legacy_reconstruction),\
         make_hlt2_tracks.bind(light_reco=False),\
-        make_reco_pvs.bind(make_pvs_from_velo_tracks=make_PatPV3DFuture_pvs),\
         make_VeloClusterTrackingSIMD.bind(algorithm=VeloClusterTrackingSIMDFull):
             #make tracks and PVs
             # TODO: in principle the PVs should be an optional input for the alignment

Alignment/Humboldt/tests/qmtest/humboldt.qms/align-ft-modules-prkalman.qmt

@@ -15,7 +15,7 @@ Verify that alignment of FT modules works with Humboldt configuration using part
 <extension class="GaudiTest.GaudiExeTest" kind="test">
   <argument name="program"><text>gaudiiter.py</text></argument>
   <argument name="args"><set>
-    <text>-n 1</text>
+    <text>-n 3</text>
     <text>-e 100</text>
     <text>$HUMBOLDTOPTS/AlignFTModules_PrKalman.py</text>
   </set></argument>

Alignment/Humboldt/tests/qmtest/humboldt.qms/align-ft-modules.qmt

@@ -15,7 +15,7 @@ Verify that alignment of FT modules works with Humboldt configuration using part
 <extension class="GaudiTest.GaudiExeTest" kind="test">
   <argument name="program"><text>gaudiiter.py</text></argument>
   <argument name="args"><set>
-    <text>-n 1</text>
+    <text>-n 3</text>
     <text>-e 100</text>
     <text>$HUMBOLDTOPTS/AlignFTModules.py</text>
   </set></argument>

Alignment/Humboldt/tests/qmtest/humboldt.qms/align-ft-stations-layers.qmt

@@ -15,7 +15,7 @@ Verify that alignment of FT stations and layers works with Humboldt configuratio
 <extension class="GaudiTest.GaudiExeTest" kind="test">
   <argument name="program"><text>gaudiiter.py</text></argument>
   <argument name="args"><set>
-    <text>-n 1</text>
+    <text>-n 3</text>
     <text>-e 100</text>
     <text>$HUMBOLDTOPTS/AlignFTStationsLayers.py</text>
   </set></argument>

Alignment/TAlignment/include/TAlignment/AlignmentElement.h

@@ -88,7 +88,7 @@ namespace LHCb::Alignment {
      */
     const std::string& fullname() const { return m_fullname; }
 
-    /** Part of the name that all detector elements have in common. May need to cash this. */
+    /** Part of the name that all detector elements have in common. */
     const std::string& basename() const { return m_basename; }
 
     /** Method to get the name of the detector(s) element
@@ -147,9 +147,6 @@ namespace LHCb::Alignment {
     /** Jacobian for the transformation from the global frame to the alignment frame */
     const Gaudi::Matrix6x6& jacobian() const { return m_jacobian; }
 
-    /** Return all elements that are served by this alignment element */
-    const ElementContainer& elementsInTree() const { return m_elementsInTree; }
-
     /** return the current delta */
     AlParameters currentTotalDelta() const;
 
@@ -261,7 +258,6 @@ namespace LHCb::Alignment {
     Gaudi::Transform3D m_alignmentFrame;                        ///< Frame in which delta-derivatives are calculated
     Gaudi::Matrix6x6   m_jacobian;                              ///< Jacobian for going from global to alignment frame
     bool               m_useLocalFrame;                         ///< Use local frame as alignmentframe
-    ElementContainer   m_elementsInTree; ///< Return all elements that are served by this alignment element
     AlParameters       m_lastDeltaDelta; ///< Last applied correction. We only keep it for its cov matrix.
   };
 

Alignment/TAlignment/include/TAlignment/GetElementsToBeAligned.h

@@ -20,6 +20,7 @@
 #ifdef USE_DD4HEP
 #  include <DD4hep/GrammarUnparsed.h>
 #endif
+#include "VPDet/DeVP.h"
 
 #include "GaudiKernel/IDataProviderSvc.h"
 #include "GaudiKernel/SmartIF.h"
@@ -72,8 +73,9 @@ namespace LHCb::Alignment {
     // Return method that finds an alignment element for a given Measuerment
     const Element* findElement( const LHCb::FitNode& node ) const;
     const Element* findElement( const LHCb::Pr::Tracks::Fit::Node& node ) const;
+    const Element* findElement( const LHCb::LHCbID id ) const;
 
-    // Find the list of elements corresponding to a path (which can ba rehulare expression)
+    // Find the list of elements corresponding to a path (which can be a regular expression)
     Elements findElements( const std::string& path ) const;
 
     // initialize the alignment frames now. if time=0, it will use the
@@ -104,26 +106,16 @@ namespace LHCb::Alignment {
     void            writeAlignmentConditions( const XmlWriter& writer, MsgStream& warningStream ) const;
     DetectorElement getDet( std::string const& path ) const;
 
-    struct ElementSorter {
-      // this is very tricky: how do you sort elements in a tree if you
-      // don't even know the tree( because of grouping). the following
-      // should cover all cases (I hope).
-
-      bool operator()( const Element* lhs, const Element* rhs ) const {
-        return ( lhs->basename() < rhs->basename() ||
-                 ( lhs->basename() == rhs->basename() &&
-                   ( lhs->detelements().front().name() < rhs->detelements().front().name() ||
-                     ( lhs->detelements().front().name() == rhs->detelements().front().name() &&
-                       lhs->elementsInTree().size() > rhs->elementsInTree().size() ) ) ) );
-      }
-    };
-
   private:
     SmartIF<IDataProviderSvc>                 m_detSvc{nullptr};
     DetectorElement                           m_topOfGeometry; /// Current top of geometry
     Elements                                  m_elements;      /// Flat vector of alignment elements
-    std::map<DetectorElement, const Element*> m_elementMap;    /// Map of detector elements to
-    bool                                      m_isInitialized{false};
+    std::map<DetectorElement, const Element*> m_elementMap;    /// Map of detector elements to alignable
+    std::map<int, const Element*>             m_vpidmap; // map from vp sensor id to alignable. FIXME: make array[0,207]
+    std::map<int, const Element*>             m_utidmap; // map from ut sensor id to alignable
+    std::map<int, const Element*>             m_ftidmap; // map from ft module id to alignable  FIXME: make array[0,257]
+    std::map<int, const Element*> m_muonidmap; // map from muon chamber id to alignable FIXME: make array[0,191]
+    bool                          m_isInitialized{false};
   };
 
 } // namespace LHCb::Alignment

Alignment/TAlignment/python/TAlignment/Alignables.py

@@ -30,7 +30,8 @@ class Alignables(list):
             "WithSupport" if not UseDD4Hep else "")
         self.m_vpModulesRight = self.m_vpSideBase + "/VPRight/Module.{1,2}" + (
             "WithSupport" if not UseDD4Hep else "")
-        self.m_vpSensors = self.m_vpSideBase + "/VP(Right|Left)/Module.{1,2}WithSupport/Module.{1,2}/Ladder."  #TODO: DD4HEP
+        self.m_vpSensors = self.m_vpModules + ("/Ladder." if not UseDD4Hep else
+                                               "/ladder_./sensor")
 
         self.m_tt = "/dd/Structure/LHCb/BeforeMagnetRegion/TT"
         self.m_ttStations = ["TTa", "TTb"]

Alignment/TAlignment/src/AlignUpdateTool.cpp

@@ -387,6 +387,9 @@ namespace LHCb::Alignment {
                << "Time at initialize [ns] : " << equations.initTime().ns() << " --> "
                << equations.initTime().format( true, "%F %r" ) << std::endl
                << "First/last run number: " << equations.firstRun() << "/" << equations.lastRun() << std::endl
+#ifdef USE_DD4HEP
+               << "Min/max IOV: " << equations.minIOV() << "/" << equations.maxIOV() << std::endl
+#endif
                << "Used " << equations.numVertices() << " vertices for alignment" << std::endl
                << "Used " << equations.numParticles() << " particles for alignment" << std::endl
                << "Used " << equations.numTracks() << " tracks for alignment" << std::endl

Alignment/TAlignment/src/AlignmentElement.cpp

@@ -24,13 +24,6 @@
 #include <iomanip>
 
 namespace {
-  void addToElementsInTree( const LHCb::Alignment::DetectorElement&     element,
-                            LHCb::Alignment::Element::ElementContainer& elements ) {
-    elements.push_back( element );
-    elements.reserve( elements.size() + element.size() );
-    element.applyToAllChildren(
-        [&]( const LHCb::Alignment::DetectorElement& child ) { addToElementsInTree( child, elements ); } );
-  }
 #ifdef USE_DD4HEP
   dd4hep::Delta amendDelta( const dd4hep::Delta& d, const Gaudi::Transform3D& m ) {
     // We should not have to do that manually ourselves, but DD4hep does not
@@ -135,9 +128,6 @@ void LHCb::Alignment::Element::addElements( const std::vector<LHCb::Alignment::D
   for ( const auto& elem : newelements )
     if ( std::find( m_elements.begin(), m_elements.end(), elem ) == m_elements.end() ) elements.push_back( elem );
 
-  for ( auto& elem : elements ) addToElementsInTree( elem, m_elementsInTree );
-  std::sort( m_elementsInTree.begin(), m_elementsInTree.end() );
-
   m_elements.insert( m_elements.end(), elements.begin(), elements.end() );
   std::sort( m_elements.begin(), m_elements.end(), SortByName() );
 
@@ -154,20 +144,6 @@ void LHCb::Alignment::Element::addElements( const std::vector<LHCb::Alignment::D
   initAlignmentFrame();
 }
 
-bool LHCb::Alignment::Element::isOffspring( const LHCb::Alignment::Element& dau ) const {
-  bool found = false;
-  if ( dau.m_elementsInTree.size() < m_elementsInTree.size() ) {
-    // check that _all_ elements are there. note that we have a real
-    // problem if only a subset is there ...
-    found = true;
-    for ( auto& elem : dau.m_elements ) {
-      found = std::find( m_elementsInTree.begin(), m_elementsInTree.end(), elem ) != m_elementsInTree.end();
-      if ( !found ) break;
-    }
-  }
-  return found;
-}
-
 void LHCb::Alignment::Element::initAlignmentFrame() {
   // Calculate the center of gravity (the 'local' center)
   double xmin( 1e10 ), xmax( -1e10 ), ymin( 1e10 ), ymax( -1e10 ), zmin( 1e10 ), zmax( -1e10 );
@@ -388,7 +364,7 @@ std::ostream& LHCb::Alignment::Element::fillStream( std::ostream& lhs ) const {
   lhs << "* Alignable: " << name() << "\n"
       << "* FullName : " << fullname() << std::endl
       << "* Element  : " << description() << std::endl
-      << "* Num elements in tree: " << m_elementsInTree.size() << std::endl
+      << "* Num elements : " << m_elements.size() << std::endl
       << "* Basename : " << basename() << std::endl
       << "* Index    : " << index() << "\n"
       << "* dPosXYZ  : " << Gaudi::XYZPoint( t[0], t[1], t[2] ) << "\n"

Alignment/TAlignment/src/GetElementsToBeAligned.cpp

@@ -126,6 +126,27 @@ namespace {
     } );
   }
 
+  /// determine if an AlignmentElement fully overlaps with det elements in a container (usually the mother)
+  template <typename DetElementContainer, typename AlignmentElement>
+  bool isOffspring( const DetElementContainer& mother_elementsInTree, const AlignmentElement& daughter ) {
+    bool found = true;
+    for ( auto& elem : daughter.detelements() ) {
+      found =
+          std::find( mother_elementsInTree.begin(), mother_elementsInTree.end(), elem ) != mother_elementsInTree.end();
+      if ( !found ) break;
+    }
+    return found;
+  }
+
+  /// create a flat container for an entire detector element tree
+  template <typename DetElementContainer>
+  void addToElementsInTree( const LHCb::Alignment::DetectorElement& element, DetElementContainer& elements ) {
+    elements.push_back( element );
+    elements.reserve( elements.size() + element.size() );
+    element.applyToAllChildren(
+        [&]( const LHCb::Alignment::DetectorElement& child ) { addToElementsInTree( child, elements ); } );
+  }
+
 } // namespace
 
 LHCb::Alignment::GetElementsToBeAligned::GetElementsToBeAligned( LHCb::span<const std::string> elemsToBeAligned,
@@ -133,8 +154,8 @@ LHCb::Alignment::GetElementsToBeAligned::GetElementsToBeAligned( LHCb::span<cons
                                                                  bool useLocalFrame, MsgStream& warningStream,
                                                                  SmartIF<IDataProviderSvc>& detSvc )
     : m_detSvc( detSvc ), m_topOfGeometry( topOfGeometry ) {
-  size_t                index( 0 );
-  std::vector<Element*> alignelements;
+  size_t index( 0 );
+  auto&  alignelements = m_elements;
   for ( const auto& elemToBeAligned : elemsToBeAligned ) {
     /// Split string into path to det elem regex and dofs to align
     bool        groupElems = false;
@@ -227,7 +248,11 @@ LHCb::Alignment::GetElementsToBeAligned::GetElementsToBeAligned( LHCb::span<cons
 
   // sort the elements by hierarchy. currently, this just follows the
   // name of the first associated detector element.
-  std::stable_sort( alignelements.begin(), alignelements.end(), ElementSorter() );
+  std::stable_sort( alignelements.begin(), alignelements.end(), []( const auto& lhs, const auto& rhs ) {
+    return ( lhs->basename() < rhs->basename() ) ||
+           ( lhs->basename() == rhs->basename() &&
+             ( lhs->detelements().front().name() < rhs->detelements().front().name() ) );
+  } );
 
   // make sure to reset the indices after sorting
   for ( auto&& [i, elem] : LHCb::range::enumerate( alignelements ) ) elem->setIndex( i );
@@ -235,10 +260,17 @@ LHCb::Alignment::GetElementsToBeAligned::GetElementsToBeAligned( LHCb::span<cons
   // set up the mother-daughter relations. This only works once the
   // elements are sorted. We start with the lowest level elements, and
   // work our way up the tree.
+
+  // first for every alignment element create a temporary flat container with all elements in the tree
+  std::vector<Alignment::Element::ElementContainer> elementsInTree( alignelements.size() );
+  for ( auto&& [i, elem] : LHCb::range::enumerate( alignelements ) )
+    for ( const auto& detelem : elem->detelements() ) addToElementsInTree( detelem, elementsInTree[i] );
+  // now use that to find the mothers
   for ( auto ielem = alignelements.rbegin(); alignelements.rend() != ielem; ++ielem ) {
     // is 'i' a daughter of 'j'
-    auto jelem = std::find_if( std::next( ielem ), alignelements.rend(),
-                               [ielem]( const auto& e ) { return e->isOffspring( **ielem ); } );
+    auto jelem = std::find_if( std::next( ielem ), alignelements.rend(), [ielem, elementsInTree]( const auto& e ) {
+      return isOffspring( elementsInTree[e->index()], **ielem );
+    } );
     if ( jelem != alignelements.rend() ) ( *jelem )->addDaughter( **ielem );
   }
 
@@ -247,11 +279,39 @@ LHCb::Alignment::GetElementsToBeAligned::GetElementsToBeAligned( LHCb::span<cons
   // alignment elements are sorted: we add the hit only to the lowest
   // level element in the tree.
   for ( const auto* ielem : alignelements ) {
-    for ( auto& detelem : ielem->elementsInTree() ) m_elementMap[detelem] = ielem;
+    for ( auto& detelem : elementsInTree[ielem->index()] ) {
+      m_elementMap[detelem] = ielem;
+#ifdef USE_DD4HEP
+      const auto* vpsensorobject =
+          dynamic_cast<const LHCb::Detector::detail::DeVPSensorObject*>( detelem.data().access() );
+      if ( vpsensorobject ) m_vpidmap[int( vpsensorobject->sensorNumber )] = ielem;
+      const auto* utsensorobject =
+          dynamic_cast<const LHCb::Detector::UT::detail::DeUTSectorObject*>( detelem.data().access() );
+      if ( utsensorobject ) m_utidmap[int( utsensorobject->m_channelID.uniqueSector() )] = ielem;
+      const auto* ftmoduleobject =
+          dynamic_cast<const LHCb::Detector::detail::DeFTModuleObject*>( detelem.data().access() );
+      if ( ftmoduleobject ) m_ftidmap[int( ftmoduleobject->m_id )] = ielem;
+      const auto* muonpadobject =
+          dynamic_cast<const LHCb::Detector::detail::DeMuonChamberObject*>( detelem.data().access() );
+      if ( muonpadobject ) m_muonidmap[int( muonpadobject->m_ChamberNumber )] = ielem;
+#else
+      const auto* vpsensorobject = dynamic_cast<const DeVPSensor*>( detelem.data() );
+      if ( vpsensorobject ) m_vpidmap[int( vpsensorobject->sensorNumber() )] = ielem;
+      const auto* ftmoduleobject = dynamic_cast<const DeFTModule*>( detelem.data() );
+      if ( ftmoduleobject ) m_ftidmap[int( ftmoduleobject->elementID().globalModuleIdx() )] = ielem;
+      const auto* utsensorobject = dynamic_cast<const DeUTSector*>( detelem.data() );
+      if ( utsensorobject ) m_utidmap[int( utsensorobject->elementID().uniqueSector() )] = ielem;
+      const auto* muonpadobject = dynamic_cast<const DeMuonChamber*>( detelem.data() );
+      if ( muonpadobject ) m_muonidmap[int( muonpadobject->chamberNumber() )] = ielem;
+#endif
+    }
   }
-
-  // copy the non-const elements into the local container.
-  m_elements.insert( m_elements.end(), alignelements.begin(), alignelements.end() );
+  // warningStream << "Elements in vpidmap:" << m_vpidmap.size() << " " << m_vpidmap.begin()->first << " " <<
+  // m_vpidmap.rbegin()->first << endmsg ; warningStream << "Elements in utidmap:" << m_utidmap.size() << " " <<
+  // m_utidmap.begin()->first << " " << m_utidmap.rbegin()->first << endmsg ; warningStream << "Elements in ftidmap:" <<
+  // m_ftidmap.size() << " " << m_ftidmap.begin()->first << " " << m_ftidmap.rbegin()->first << endmsg ; warningStream
+  // << "Elements in muonidmap:" << m_muonidmap.size() << " " << m_muonidmap.begin()->first << " " <<
+  // m_muonidmap.rbegin()->first << endmsg ;
 }
 
 LHCb::Alignment::GetElementsToBeAligned::~GetElementsToBeAligned() {
@@ -299,28 +359,42 @@ LHCb::Alignment::Elements LHCb::Alignment::GetElementsToBeAligned::findElements(
 
 const LHCb::Alignment::Element*
 LHCb::Alignment::GetElementsToBeAligned::findElement( const LHCb::FitNode& node ) const {
-  const Element* elem = nullptr;
-  try {
-    if ( node.hasMeasurement() ) {
-      auto detelem = node.measurement().visit(
-#ifdef USE_DD4HEP
-          []( const auto& ) -> DetectorElement { throw std::logic_error( "Subdetector not supported in DD4hep" ); },
-#else
-          []( const auto& m ) { return DetectorElement( m.detElem ); },
-#endif
-          []( const Measurement::VP& vp ) { return DetectorElement( vp.detElem ); } );
-      return findElement( detelem );
-    }
-  } catch ( std::logic_error& ) {
-    // FIXME Ignoring parts not supported in DD4hep
-  }
-  return elem;
+  return findElement( node.measurement().lhcbID() );
 }
 
 const LHCb::Alignment::Element*
 LHCb::Alignment::GetElementsToBeAligned::findElement( const LHCb::Pr::Tracks::Fit::Node& node ) const {
-  auto detelem = m_topOfGeometry.findElemForPoint( {node.final_state_vec[0], node.final_state_vec[1], node.z()} );
-  return detelem ? findElement( *detelem ) : nullptr;
+  return findElement( node.lhcbID );
+}
+
+const LHCb::Alignment::Element* LHCb::Alignment::GetElementsToBeAligned::findElement( const LHCb::LHCbID id ) const {
+  const Element* element{0};
+  switch ( id.detectorType() ) {
+  case LHCb::LHCbID::channelIDtype::VP: {
+    auto it = m_vpidmap.find( int( id.vpID().sensor() ) );
+    if ( it != m_vpidmap.end() ) element = it->second;
+  } break;
+  case LHCb::LHCbID::channelIDtype::UT: {
+    // FIXME: At some point we need here also the element of the trajectory: see run2_patches
+    auto it = m_utidmap.find( int( id.utID().uniqueSector() ) );
+    if ( it != m_utidmap.end() ) element = it->second;
+  } break;
+  case LHCb::LHCbID::channelIDtype::FT: {
+    auto it = m_ftidmap.find( int( id.ftID().globalModuleIdx() ) );
+    if ( it != m_ftidmap.end() ) element = it->second;
+  } break;
+  case LHCb::LHCbID::channelIDtype::Muon: {
+    auto it = m_muonidmap.find( int( id.muonID() ) );
+    if ( it != m_muonidmap.end() ) element = it->second;
+  } break;
+    // case LHCb::LHCbID::Muon:
+    //   element = muon->getChmbPtr( id.muonID().station(), id.muonID().region(),
+    // 				( muon->Tile2Chamber( id.muonID() ) )[0]->chamberNumber() );
+    //   break;
+  default: {
+  };
+  }
+  return element;
 }
 
 const LHCb::Alignment::Element*

Alignment/TAlignment/src/TrackResidualHelper.cpp

@@ -302,12 +302,12 @@ namespace LHCb::Alignment {
         residuals->m_HCHElements.reserve( numalignables * ( numalignables - 1 ) / 2 );
         for ( size_t irow = 0; irow < numalignables; ++irow ) {
           size_t k = alignable_node_indices[irow];
-          assert( &( residuals->m_residuals[irow].node() ) == nodeswithindex[k].first );
+          assert( std::get<const TNode*>( residuals->m_residuals[irow].node() ) == nodeswithindex[k].first );
           const Gaudi::TrackProjectionMatrix1D& Hk = projectionMatrix( nodeswithindex[k].first );
           // first the off-diagonal elements
           for ( size_t icol = 0; icol < irow; ++icol ) {
             size_t l = alignable_node_indices[icol];
-            assert( &( residuals->m_residuals[icol].node() ) == nodeswithindex[l].first );
+            assert( std::get<const TNode*>( residuals->m_residuals[icol].node() ) == nodeswithindex[l].first );
             const Gaudi::TrackProjectionMatrix1D& Hl  = projectionMatrix( nodeswithindex[l].first );
             double                                hch = l < k
                              ? ( Hk * getOffDiagCov( k, l, offdiagcov, smoothergainmatrix ) * Transpose( Hl ) )( 0, 0 )