Adapt to updated UT channel ID class and geometry

(1) Changes based on new UT geometry & channelID (2) stripID shifted. Now start from 0 instead of 1

Work with DDDB!79, SIMCOND!143, LHCb!2995 (merged), Lbcom!550 (merged), Moore!785 (merged), Allen!865 (merged), PRConfig!230, Detector!181 (merged)

Pr/PrMCTools/src/PrUTHitsMonitor.cpp

+/*****************************************************************************\
+ * * (c) Copyright 2000-2020 CERN for the benefit of the LHCb Collaboration      *
+ * *                                                                             *
+ * * This software is distributed under the terms of the GNU General Public      *
+ * * Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
+ * *                                                                             *
+ * * In applying this licence, CERN does not waive the privileges and immunities *
+ * * granted to it by virtue of its status as an Intergovernmental Organization  *
+ * * or submit itself to any jurisdiction.                                       *
+ * \*****************************************************************************/
+#include "DetDesc/GenericConditionAccessorHolder.h"
+#include "DetDesc/IConditionDerivationMgr.h"
+#include "Event/MCHit.h"
+#include "GaudiAlg/Consumer.h"
+#include "GaudiAlg/GaudiHistoAlg.h"
+#include "GaudiKernel/Vector3DTypes.h"
+#include "Kernel/LHCbID.h"
+#include "UTDAQ/UTInfo.h"
+
+#include "AIDA/IHistogram2D.h"
+#include "AIDA/IProfile1D.h"
+#include "Event/LinksByKey.h"
+#include "Event/MCParticle.h"
+#include "Event/UTCluster.h"
+#include "Event/UTDigit.h"
+#include "GaudiUtils/Aida2ROOT.h"
+#include "PrKernel/PrUTHitHandler.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "UTDet/DeUTDetector.h"
+
+//-----------------------------------------------------------------------------
+// Implementation file for class : PrUTHitsMonitor
+//
+// 2021-05: valeriia lukashenko
+//-----------------------------------------------------------------------------
+// @class PrUTHitsMonitor PrUTHitsMonitor.h
+// 	    * class for plotting redisuals from UTHits to MCHits
+//
+// 	    * isGlobal - if true use global coordinates, if false use local coordinates
+// 	    * performStudy  - if true plot residuals, residulas in different planes,
+// residuals vs track angle. 	   		      if false just plot residuals
+// 	    * min is the left side of the residuals range
+// 	    * max is the right side of the residuals range
+// 	    * nbins is number of the bins for residuals
+//          * to be used in the Dashboard
+//          * example: Moore/Hlt/RecoCond/options/tracking_developments/mc_hit_resolution_monitor.py
+//
+//
+// *  @author valeriia lukashenko
+// *  @date   2021-05-04
+//
+namespace {
+  using MCHits      = LHCb::MCHits;
+  using MCHit       = LHCb::MCHit;
+  using MCParticles = LHCb::MCParticles;
+  using UTHits      = LHCb::Pr::UT::HitHandler;
+  using simd        = SIMDWrapper::scalar::types;
+} // namespace
+class PrUTHitsMonitor final
+    : public Gaudi::Functional::Consumer<void( const UTHits&, const LHCb::LinksByKey&, const MCHits&,
+                                               const MCParticles&, const DeUTDetector& ),
+                                         LHCb::DetDesc::usesBaseAndConditions<GaudiHistoAlg, DeUTDetector>> {
+public:
+  PrUTHitsMonitor( const std::string& name, ISvcLocator* pSvcLocator )
+      : Consumer( name, pSvcLocator,
+                  {KeyValue{"UTHitsLocation", UT::Info::HitLocation},
+                   KeyValue{"UTHits2MCHitLinksLocation", LHCb::UTDigitLocation::UTTightDigits + "2MCHits"},
+                   //KeyValue{"UTHits2MCHitLinksLocation", LHCb::UTClusterLocation::UTClusters + "2MCHits"},
+                   KeyValue{"MCHitsLocation", LHCb::MCHitLocation::UT},
+                   KeyValue{"MCParticleLocation", LHCb::MCParticleLocation::Default},
+                   KeyValue{"DeUT", DeUTDetLocation::UT}} ) {}
+
+  void                    operator()( const UTHits&, const LHCb::LinksByKey&, const MCHits&, const MCParticles&,
+                   const DeUTDetector& ) const override;
+  Gaudi::Property<bool>   m_isGlobal{this, "isGlobal", true, "Boolean for choosing global coordinate system"};
+  Gaudi::Property<bool>   m_performStudy{this, "performStudy", false,
+                                       "Boolean for performing studies with cluster size and diraction of mcparticle"};
+  Gaudi::Property<double> m_min{this, "m_min_bin", -0.5, "Double for left boundary of residuals range"};
+  Gaudi::Property<double> m_max{this, "m_max_bin", 0.5, "Double for right boundary of residuals range"};
+  Gaudi::Property<int>    m_num_bins{this, "m_num_bins", 1000, "Int for number of bins"};
+};
+
+DECLARE_COMPONENT( PrUTHitsMonitor )
+
+void PrUTHitsMonitor::operator()( const UTHits& UThits, const LHCb::LinksByKey& links, const MCHits& mchits,
+                                  const MCParticles&, const DeUTDetector& DeUT ) const {
+
+  std::map<const unsigned int, std::vector<LHCb::MCHit*>> mcHitForId;
+  links.applyToAllLinks( [&mcHitForId, &mchits]( unsigned int id, unsigned int mcHitKey, float ) {
+    mcHitForId[id].emplace_back( mchits[mcHitKey] );
+  } );
+
+#if 0
+  //To Be Done
+  // Map for MCHits and keys
+  if ( !m_performStudy ) {
+    const auto& hits = UThits.hits();
+    const int   fullChanIdx =
+        static_cast<int>( UTInfo::DetectorNumbers::Layers ) * static_cast<int>( UTInfo::DetectorNumbers::Stations ) *
+        static_cast<int>( UTInfo::DetectorNumbers::Regions ) * static_cast<int>( UTInfo::DetectorNumbers::Sectors );
+    // new SOA for UTHits: hit.get<LHCb::Pr::UT::UTHitsTag::zAtYEq0>
+    for ( int fullchan = 0; fullchan < fullChanIdx; fullchan++ ) {
+      const auto indexs = UThits.indices( fullchan );
+
+      for ( int i = indexs.first; i != indexs.second; i++ ) {
+        const auto        hit           = hits.scalar()[i];
+        const simd::int_v simd_chid     = hit.get<LHCb::Pr::UT::UTHitsTag::channelID>();
+        const int         unwrapped_num = simd_chid.cast();
+
+        for ( auto& mcHit : ( *mcHitForId.find( unwrapped_num ) ).second ) {
+          const auto hit_x = hit.get<LHCb::Pr::UT::UTHitsTag::xAtYEq0>().cast() +
+                             mcHit->midPoint().Y() * hit.get<LHCb::Pr::UT::UTHitsTag::dxDy>().cast();
+          Gaudi::XYZPoint UThit = Gaudi::XYZPoint( hit_x, 0, hit.get<LHCb::Pr::UT::UTHitsTag::zAtYEq0>().cast() );
+
+          auto residual = UThit - ( mcHit->midPoint() );
+          plot1D( residual.X(), "residuals_X", m_min, m_max, m_num_bins );
+          plot1D( residual.Z(), "residuals_Z", m_min, m_max, m_num_bins );
+
+        } // while mcHit
+      }   // hit loop
+    }     // channel id loop*/
+  }       // if not  performStudy
+  if ( m_performStudy ) {
+    AIDA::IHistogram2D* res_vs_p_x       = nullptr;
+    AIDA::IProfile1D*   tprof_rms_vs_p_x = nullptr;
+    if ( m_isGlobal ) {
+      res_vs_p_x = book2D( "residuals_X_VS_slope_dxdz", m_min, m_max, m_num_bins / 10., -1., 1., m_num_bins / 10. );
+      tprof_rms_vs_p_x = bookProfile1D( "tprof_res_rms_vs_slope_dxdz", -1., 1., 100. );
+    } else {
+      res_vs_p_x =
+          book2D( "residuals_local_X_VS_slope_dxdz", m_min, m_max, m_num_bins / 10., -1., 1., m_num_bins / 10. );
+      tprof_rms_vs_p_x = bookProfile1D( "tprof_res_rms_local_vs_slope_dxdz", -1., 1., 100. );
+    }
+
+    const auto& hits = UThits.hits();
+    const int   fullChanIdx =
+        static_cast<int>( UTInfo::DetectorNumbers::Layers ) * static_cast<int>( UTInfo::DetectorNumbers::Stations ) *
+        static_cast<int>( UTInfo::DetectorNumbers::Regions ) * static_cast<int>( UTInfo::DetectorNumbers::Sectors );
+
+    for ( int fullchan = 0; fullchan < fullChanIdx; fullchan++ ) {
+      const auto indexs = UThits.indices( fullchan );
+
+      for ( int i = indexs.first; i != indexs.second; i++ ) {
+        const auto        hit           = hits.scalar()[i];
+        const simd::int_v simd_chid     = hit.get<LHCb::Pr::UT::UTHitsTag::channelID>().cast();
+        const int         unwrapped_num = simd_chid.cast();
+        auto              sector        = DeUT.findSector( LHCb::UTChannelID( unwrapped_num ) );
+        auto              result        = ( mcHitForId.find( unwrapped_num ) );
+        if ( result == mcHitForId.end() ) continue;
+        for ( auto& mcHit : ( *result ).second ) {
+          const auto hit_x = hit.get<LHCb::Pr::UT::UTHitsTag::xAtYEq0>().cast() +
+                             mcHit->midPoint().Y() * hit.get<LHCb::Pr::UT::UTHitsTag::dxDy>().cast();
+          Gaudi::XYZPoint   UThit = Gaudi::XYZPoint( hit_x, 0, hit.get<LHCb::Pr::UT::UTHitsTag::zAtYEq0>().cast() );
+          LHCb::UTChannelID utid  = LHCb::UTChannelID( (unsigned int)unwrapped_num );
+          if ( !sector )
+            error() << "No sector is found for "
+                    << "(" << hit.get<LHCb::Pr::UT::UTHitsTag::xAtYEq0>().cast() << ", 0, "
+                    << hit.get<LHCb::Pr::UT::UTHitsTag::zAtYEq0>().cast() << ")" << std::endl;
+          const auto UThit_local = sector->toLocal( UThit );
+          const auto mcHit_local = sector->toLocal( mcHit->midPoint() );
+
+          const int        plane = 2 * ( utid.station() - 1 ) + ( utid.layer() - 1 ) % 2;
+          Gaudi::XYZVector residual;
+          std::string      name_addition = "";
+          if ( m_isGlobal ) {
+            residual = UThit - ( mcHit->midPoint() );
+          } else {
+            residual      = UThit_local - mcHit_local;
+            name_addition = "_local";
+          }
+
+          plot1D( residual.X(), "residuals" + name_addition + "_X", m_min, m_max, m_num_bins );
+          plot1D( residual.Z(), "residuals" + name_addition + "_Z", m_min, m_max, m_num_bins );
+
+          plot1D( residual.X(), "residuals" + name_addition + "_X_plane" + std::to_string( ( plane ) ), m_min * 3,
+                  m_max * 3, m_num_bins );
+          plot1D( residual.Z(), "residuals" + name_addition + "_Z_plane" + std::to_string( ( plane ) ), m_min, m_max,
+                  m_num_bins );
+          auto       MCparticle = mcHit->mcParticle();
+          const auto momentum   = MCparticle->momentum();
+          if ( momentum.Z() != 0 ) {
+            plot2D( residual.X(), ( momentum.X() / momentum.Z() ), "residuals" + name_addition + "_X_vs_track_angle",
+                    m_min, m_max, -1, 1, m_num_bins, m_num_bins / 10. );
+            profile1D( residual.X(), momentum.X() / momentum.Z(), "profile_res" + name_addition + "_X_VS_slope_dxdz",
+                       m_min, m_max, m_num_bins / 10. );
+            res_vs_p_x->fill( residual.X(), momentum.X() / momentum.Z() );
+
+          } else {
+            plot2D( residual.X(), -9999., "residuals" + name_addition + "_X_vs_track_angle", m_min, m_max, -1, 1,
+                    m_num_bins, m_num_bins / 10. );
+            profile1D( residual.X(), -9999., "profile_res" + name_addition + "_X_VS_slope_dxdz", m_min, m_max,
+                       m_num_bins / 10. );
+
+            res_vs_p_x->fill( residual.X(), -9999. );
+          }
+          plot1D( residual.X(), "residuals" + name_addition + "_X_plane_" + std::to_string( plane ) + "_zoom", m_min,
+                  m_max, m_num_bins );
+        } // loop mcHits
+      }   // loop indices
+      TH2D* th2d_res_vs_p_x = Gaudi::Utils::Aida2ROOT::aida2root( res_vs_p_x );
+      if ( !m_isGlobal ) { th2d_res_vs_p_x->SetName( "residuals_local_X_VS_track_slope_momentum" ); }
+      for ( int bin = 0; bin != th2d_res_vs_p_x->GetNbinsX(); bin++ ) {
+        TH1D*  th1d_x    = th2d_res_vs_p_x->ProjectionX( "th1d_x", bin, bin );
+        double RMS_x     = th1d_x->GetRMS();
+        double bin_value = th1d_x->GetXaxis()->GetBinCenter( bin );
+        tprof_rms_vs_p_x->fill( bin_value, RMS_x );
+      }
+
+    } // if performStudy
+  }
+#endif
+}