MuonSpacepointformation - Fix covariance & misc

MuonSpacepointformation - Fix covariance & misc

MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/MuonSpacePoint.cxx

@@ -53,7 +53,8 @@ namespace MuonR4{
             }
             uvcov(1,1) = std::pow(uvcov(1,1), 2);
         }
-        AmgSymMatrix(2) cov = Jac.inverse() * uvcov * Jac; 
+        Jac = Jac.inverse();
+        AmgSymMatrix(2) cov = Jac * uvcov * Jac.transpose();
         m_measUncerts =  Amg::Vector2D(std::sqrt(cov(0,0)), std::sqrt(cov(1,1)));
     }
             

MuonSpectrometer/MuonPhaseII/Event/xAOD/xAODMuonPrepData/Root/RpcStrip_v1.cxx

@@ -47,4 +47,3 @@ ConstVectorMap<3> RpcStrip_v1::stripPosInStation() const {
 
 }  // namespace xAOD
 #undef IMPLEMENT_SETTER_GETTER
-#undef THROW_EXCEPT

MuonSpectrometer/MuonPhaseII/Event/xAOD/xAODMuonPrepData/Root/UtilFunctions.cxx

@@ -2,7 +2,7 @@
    Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 #include "xAODMuonPrepData/UtilFunctions.h"
-
+#include "GeoModelHelpers/throwExcept.h"
 #include "xAODMuonPrepData/MdtDriftCircle.h"
 #include "xAODMuonPrepData/RpcStrip.h"
 #include "xAODMuonPrepData/TgcStrip.h"
@@ -40,7 +40,7 @@ namespace xAOD{
         } else if (meas->type() == xAOD::UncalibMeasType::sTgcStripType) {
             return static_cast<const xAOD::sTgcMeasurement*>(meas)->readoutElement();
         }
-        THROW_EXCEPT("Unsupported measurement given "<<typeid(*meas).name());
+        THROW_EXCEPTION("Unsupported measurement given "<<typeid(*meas).name());
         return nullptr;
     }
     
@@ -59,12 +59,12 @@ namespace xAOD{
             const xAOD::TgcStrip* strip = static_cast<const xAOD::TgcStrip*>(meas);
             return toChamberTransform(gctx, strip) *(strip->localPosition<1>()[0] * Amg::Vector3D::UnitX());
         } else {
-            THROW_EXCEPT("Measurement "<<typeid(*meas).name()<<" is not supported");
+            THROW_EXCEPTION("Measurement "<<typeid(*meas).name()<<" is not supported");
         }
         return Amg::Vector3D::Zero();
     }
     Amg::Vector3D channelDirInChamber(const ActsGeometryContext& gctx,
-                                    const xAOD::UncalibratedMeasurement* meas) {        
+                                      const xAOD::UncalibratedMeasurement* meas) {        
         if (!meas) return Amg::Vector3D::Zero();
         if (meas->type() == xAOD::UncalibMeasType::MdtDriftCircleType) {
             const xAOD::MdtDriftCircle* dc = static_cast<const xAOD::MdtDriftCircle*>(meas);
@@ -81,7 +81,7 @@ namespace xAOD{
             } 
             return trf.linear() *dir;
         }
-        THROW_EXCEPT("Measurement "<<typeid(*meas).name()<<" is not supported");
+        THROW_EXCEPTION("Measurement "<<typeid(*meas).name()<<" is not supported");
         return Amg::Vector3D::Zero();        
     }
     Amg::Vector3D channelNormalInChamber(const ActsGeometryContext& gctx,
@@ -102,7 +102,7 @@ namespace xAOD{
             } 
             return trf.linear() *dir;
         }
-        THROW_EXCEPT("Measurement "<<typeid(*meas).name()<<" is not supported");
+        THROW_EXCEPTION("Measurement "<<typeid(*meas).name()<<" is not supported");
         return Amg::Vector3D::Zero();  
     }
 }

MuonSpectrometer/MuonPhaseII/Event/xAOD/xAODMuonPrepData/xAODMuonPrepData/versions/AccessorMacros.h

@@ -3,6 +3,7 @@
 */
 #ifndef XAODMUONPREPDATA_ACCESSOR_MACROS_H
 #define XAODMUONPREPDATA_ACCESSOR_MACROS_H
+#include "GeoModelHelpers/throwExcept.h"
 /**
  *  Macros to implement the scalar variables of the xAOD::MuonPrepData objects
 */
@@ -43,17 +44,6 @@
          static const SG::AuxElement::Accessor<std::vector<DTYPE>> acc{preFixStr + #GETTER}; \
          acc(*this) = value;                                                                 \
       }
-/**
- *  Macro to throw an exception message indicating the
- *  location in the source code where the exception has been thrown
-*/
-#define THROW_EXCEPT(MSG)                       \
-    {                                           \
-       std::stringstream sstr{};                \
-       sstr<<__FILE__<<":"<<__LINE__<<" "<<MSG; \
-       throw std::runtime_error(sstr.str());    \
-    } 
-
 /**
  *  Macro to handle the readoutElement. If the object is created within the RDO -> Prd conversion
  *  the method simply returns the pointer to the given readoutElement. In contrast, if the object is 
@@ -73,14 +63,14 @@
             const MuonGMR4::MuonDetectorManager* detMgr{};                                      \
             if (!service.retrieve().isSuccess() ||                                              \
                 !service->retrieve(detMgr).isSuccess()){                                        \
-                THROW_EXCEPT("Failed to retrieve the Run4 muon detector manager. "<<            \
+                THROW_EXCEPTION("Failed to retrieve the Run4 muon detector manager. "<<         \
                          "Please schedule the MuonGeometry in your job");                       \
             }                                                                                   \
             const IdentifierHash hash{identifierHash()};                                        \
             const MuonGMR4::READOUT_ELEMENT_TYPE* re = detMgr->get##READOUT_ELEMENT_TYPE(hash); \
             if (!re) {                                                                          \
                 const Identifier id{static_cast<Identifier::value_type>(identifier())};         \
-                THROW_EXCEPT(detMgr->idHelperSvc()->toString(id)                                \
+                THROW_EXCEPTION(detMgr->idHelperSvc()->toString(id)                             \
                             <<" does not have a readout element.");                             \
             }                                                                                   \
             CACHED_VALUE.set(re);                                                               \

MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughHelpers.cxx

@@ -5,27 +5,27 @@
 #include "MuonHoughHelpers.h"
 using namespace MuonR4; 
 
-double HoughHelpers::Eta::houghParamMdtLeft(double tanTheta, const MuonR4::HoughHitType & DC){
+  double HoughHelpers::Eta::houghParamMdtLeft(double tanTheta, const MuonR4::HoughHitType & DC){
     return DC->positionInChamber().y() - tanTheta * DC->positionInChamber().z() -
            DC->driftRadius() * std::sqrt(1 + (tanTheta*tanTheta));    // using cos(theta) = sqrt(1/[1+tan²(theta)])
-  };
+  }
   // right solution
-double HoughHelpers::Eta::houghParamMdtRight(double tanTheta, const MuonR4::HoughHitType & DC){
-  return DC->positionInChamber().y() - tanTheta * DC->positionInChamber().z() +
+  double HoughHelpers::Eta::houghParamMdtRight(double tanTheta, const MuonR4::HoughHitType & DC){
+    return DC->positionInChamber().y() - tanTheta * DC->positionInChamber().z() +
            DC->driftRadius() * std::sqrt(1 + (tanTheta*tanTheta));    // using cos(theta) = sqrt(1/[1+tan²(theta)])
-  };
+  }
   // solution which doesn't look at the drift circle but places a hit at the tube center
-double HoughHelpers::Eta::houghParamStrip(double tanTheta, const MuonR4::HoughHitType & DC){
+  double HoughHelpers::Eta::houghParamStrip(double tanTheta, const MuonR4::HoughHitType & DC){
     return DC->positionInChamber().y() - tanTheta * DC->positionInChamber().z();
-  };
+  }
 
-double HoughHelpers::Eta::houghWidthMdt(double /*tanTheta*/, const MuonR4::HoughHitType & DC){
+  double HoughHelpers::Eta::houghWidthMdt(double /*tanTheta*/, const MuonR4::HoughHitType & DC){
     return std::min(DC->uncertainty().x() * 3.,
                     1.0);  // scale reported errors up to at least 1mm or 3
                            // times the reported error as drift circle calib not
                            // fully reliable at this stage
-  };
-double HoughHelpers::Eta::houghWidthStrip(double /*tanTheta*/, const MuonR4::HoughHitType & DC){
-    return DC->uncertainty().x();  // assign full tube radius as uncertainty
-  };
+  }
+  double HoughHelpers::Eta::houghWidthStrip(double /*tanTheta*/, const MuonR4::HoughHitType & DC){
+      return DC->uncertainty().y();  // assign full tube radius as uncertainty
+  }
 

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/python/MdtEtaTransformTesterConfig.py

@@ -14,7 +14,14 @@ if __name__=="__main__":
     from MuonGeoModelTestR4.testGeoModel import setupGeoR4TestCfg, SetupArgParser, executeTest,setupHistSvcCfg
     parser = SetupArgParser()
     parser.set_defaults(nEvents = -1)
+    parser.set_defaults(noMM=True)
+    parser.set_defaults(noSTGC=True)
     parser.set_defaults(inputFile=["/cvmfs/atlas-nightlies.cern.ch/repo/data/data-art/MuonRecRTT/R4SimHits.pool.root"])
+    parser.add_argument("--displayFailedSeeds", 
+                        help="Saves the hits of failed seeds in a pdf", action='store_true', default = False)
+    parser.add_argument("--displayGoodSeeds", 
+                        help="Saves the hits of failed seeds in a pdf", action='store_true', default = False)
+
 
     args = parser.parse_args()
     flags, cfg = setupGeoR4TestCfg(args)
@@ -28,7 +35,9 @@ if __name__=="__main__":
     from MuonSpacePointFormation.SpacePointFormationConfig import MuonSpacePointMakerAlgCfg 
     cfg.merge(MuonSpacePointMakerAlgCfg(flags))
     cfg.merge(MuonHoughTransformAlgCfg(flags))
-    cfg.merge(MdtEtaTransformTesterCfg(flags))
+    cfg.merge(MdtEtaTransformTesterCfg(flags,
+                                       drawDisplayFailed =args.displayFailedSeeds,
+                                       drawDisplaySuccss = args.displayGoodSeeds))
     cfg.merge(PerfMonMTSvcCfg(flags))
     # cfg.merge(VTuneProfilerServiceCfg(flags, ProfiledAlgs=["MuonHoughTransformAlg"]))
 

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MdtEtaTransformTester.cxx

@@ -3,6 +3,7 @@
 */
 
 #include "MdtEtaTransformTester.h"
+#include "GaudiKernel/SystemOfUnits.h"
 #include "MuonReadoutGeometryR4/MuonChamber.h"
 #include "StoreGate/ReadCondHandle.h"
 #include "GeoModelHelpers/throwExcept.h"
@@ -33,11 +34,16 @@ namespace MuonValR4 {
         ATH_CHECK(m_idHelperSvc.retrieve());
         ATH_CHECK(detStore()->retrieve(m_r4DetMgr));
         ATH_MSG_DEBUG("Succesfully initialised");
+        if (m_drawEvtDisplayFailure || m_drawEvtDisplaySuccess) {
+            m_allCan = std::make_unique<TCanvas>("AllDisplays", "AllDisplays", 800, 600);
+            m_allCan->SaveAs(Form("%s[", m_allCanName.value().c_str()));
+        }
         return StatusCode::SUCCESS;
     }
 
     StatusCode MdtEtaTransformTester::finalize() {
         ATH_CHECK(m_tree.write());
+        if (m_allCan) m_allCan->SaveAs(Form("%s]", m_allCanName.value().c_str()));
         return StatusCode::SUCCESS;
     }
 
@@ -54,15 +60,15 @@ namespace MuonValR4 {
 
     StatusCode MdtEtaTransformTester::execute()  {
         
-        const EventContext & context = Gaudi::Hive::currentContext();
-        SG::ReadCondHandle<ActsGeometryContext> gctxHandle{m_geoCtxKey, context};
+        const EventContext & ctx = Gaudi::Hive::currentContext();
+        SG::ReadCondHandle<ActsGeometryContext> gctxHandle{m_geoCtxKey, ctx};
         ATH_CHECK(gctxHandle.isValid());
         const ActsGeometryContext& gctx{**gctxHandle};
         // retrieve the two input collections
 
-        auto simHitCollections = m_inSimHitKeys.makeHandles(context);
+        auto simHitCollections = m_inSimHitKeys.makeHandles(ctx);
         
-        SG::ReadHandle<MuonR4::StationHoughMaxContainer> readHoughPeaks(m_inHoughMaximaKey, context);
+        SG::ReadHandle<MuonR4::StationHoughMaxContainer> readHoughPeaks(m_inHoughMaximaKey, ctx);
         ATH_CHECK(readHoughPeaks.isPresent());        
 
         ATH_MSG_DEBUG("Succesfully retrieved input collections");
@@ -105,7 +111,7 @@ namespace MuonValR4 {
             const std::optional<double> lambda = Amg::intersect<3>(localPos, chamberDir, Amg::Vector3D::UnitZ(), 0.);
             Amg::Vector3D chamberPos = localPos + (*lambda)*chamberDir;
 
-            m_evtNumber = context.eventID().event_number();
+            m_evtNumber = ctx.eventID().event_number();
             m_out_stationName = reElement->stationName();
             m_out_stationEta = reElement->stationEta();
             m_out_stationPhi = reElement->stationPhi();
@@ -129,7 +135,12 @@ namespace MuonValR4 {
 
             const std::vector<MuonR4::HoughMaximum>& houghMaxima = houghPeakMap[reElement->getChamber()];
             if (houghMaxima.empty()){
-                if (!m_tree.fill(context)) return StatusCode::FAILURE; 
+                if (!m_tree.fill(ctx)) {
+                    return StatusCode::FAILURE;
+                }
+                if (m_drawEvtDisplayFailure) {
+                    ATH_CHECK(drawEventDisplay(ctx, hits, nullptr));
+                }
                 continue;
             }
             const MuonR4::HoughMaximum* foundMax = nullptr; 
@@ -202,9 +213,13 @@ namespace MuonValR4 {
                 m_out_max_nMdt = nMdt;
                 m_out_max_nRpc = nRpc;
                 m_out_max_nTgc = nTgc;
-            }           
- 
-            if (!m_tree.fill(context)) return StatusCode::FAILURE;
+                if (m_drawEvtDisplaySuccess) {
+                    ATH_CHECK(drawEventDisplay(ctx, hits, foundMax));
+                }
+            } else if (m_drawEvtDisplayFailure) {
+                ATH_CHECK(drawEventDisplay(ctx, hits, nullptr));          
+            }
+            if (!m_tree.fill(ctx)) return StatusCode::FAILURE;
         }
         return StatusCode::SUCCESS;
     }
@@ -277,7 +292,9 @@ namespace MuonValR4 {
                 ATH_MSG_DEBUG( "         HIT @ "<<Amg::toString(hit->positionInChamber())<<"  "<<m_idHelperSvc->toString(hit->identify())<<" with r = "<<hit->driftRadius()); 
                 /// Space point is a mdt space point
                 if (hit->driftRadius() > 1e-6) {
-                    auto  ell = std::make_unique<TEllipse>(hit->positionInChamber().y(), hit->positionInChamber().z(),hit->driftRadius());
+                    auto  ell = std::make_unique<TEllipse>(hit->positionInChamber().y(), 
+                                                          hit->positionInChamber().z(),
+                                                          hit->driftRadius());
                     ell->SetLineColor(kRed);
                     ell->SetFillColor(kRed);
                     if (hitsOnMax.count(hit)) {
@@ -289,7 +306,9 @@ namespace MuonValR4 {
                     ell->Draw(); 
                     primitives.emplace_back(std::move(ell));
                 } else {
-                    auto  m = std::make_unique<TEllipse>(hit->positionInChamber().y(), hit->positionInChamber().z(), hit->uncertainty().x(), hit->uncertainty().x());
+                    auto  m = std::make_unique<TEllipse>(hit->positionInChamber().y(), 
+                                                         hit->positionInChamber().z(), 
+                                                         hit->uncertainty().y());
                     m->SetLineColor(kBlue); 
                     m->SetFillColor(kBlue); 
                     m->SetFillStyle(0);
@@ -300,7 +319,7 @@ namespace MuonValR4 {
                     } else if (hit->measuresPhi())  {
                         m->SetLineColor(kGreen); 
                         m->SetFillColor(kGreen);
-                    } 
+                    }
                     /// Fill the ellipse if the hit is on the hough maximum
                     if (hitsOnMax.count(hit)) {
                         m->SetFillStyle(1001);
@@ -341,6 +360,7 @@ namespace MuonValR4 {
                 <<refChamber->stationEta()<<"_"<<refChamber->stationPhi()<<".pdf"; 
         
         myCanvas.SaveAs(pdfName.str().c_str());
+        myCanvas.SaveAs(m_allCanName.value().c_str());
         primitives.clear();
         return StatusCode::SUCCESS;
     }

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MdtEtaTransformTester.h

@@ -25,7 +25,7 @@
 #include "MuonTesterTree/ThreeVectorBranch.h"
 #include "MuonTesterTree/IdentifierBranch.h"
 
-
+#include "TCanvas.h"
 ///  @brief Lightweight algorithm to read xAOD MDT sim hits and 
 ///  (fast-digitised) drift circles from SG and fill a 
 ///  validation NTuple with identifier and drift circle info.
@@ -104,6 +104,14 @@ namespace MuonValR4{
     MuonVal::VectorBranch<float>&   m_max_tgcHitErrorX{m_tree.newVector<float>("maxTgcEtaMeasError")};
     MuonVal::VectorBranch<float>&   m_max_tgcHitErrorY{m_tree.newVector<float>("maxTgcPhiMeasError")};
 
+    /// Draw the event display for the cases where the hough transform did not find any hough maximum
+    Gaudi::Property<bool> m_drawEvtDisplayFailure{this, "drawDisplayFailed", false};
+    /// Draw the event dispalty for the successful cases
+    Gaudi::Property<bool> m_drawEvtDisplaySuccess{this, "drawDisplaySuccss", false};
+    
+    std::unique_ptr<TCanvas> m_allCan{};
+    Gaudi::Property<std::string> m_allCanName{this, "AllCanvasName", "AllHoughiDiPuffDisplays.pdf"};
+
   };
 }