Add a validation package for the Phase-2 muon fast MDT digitisation

Add a validation package for the Phase-2 muon fast MDT digitisation

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonFastDigiTest/CMakeLists.txt

+################################################################################
+# Package: MuonFastDigiTest
+################################################################################
+
+# Declare the package name:
+atlas_subdir( MuonFastDigiTest )
+
+atlas_add_component( MuonFastDigiTest
+                     src/components/*.cxx src/*.cxx
+                     LINK_LIBRARIES AthenaKernel StoreGateLib MuonTesterTreeLib
+                                    xAODMuonSimHit xAODMuonPrepData  )
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonFastDigiTest/python/MuonFastDigiTesterConfig.py

+# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+
+def MDTFastDigiTesterCfg(flags, name = "MDTFastDigiTester", **kwargs):
+    result = ComponentAccumulator()
+    theAlg = CompFactory.MuonValR4.MDTFastDigiTester(name, **kwargs)
+    result.addEventAlgo(theAlg, primary=True)
+    return result

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonFastDigiTest/python/testFastMdtDigi.py

+# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+
+# Small CA snippet to run the MDT fast-digi test on an xAOD file containing sim hits.
+# Will schedule fast digitisation and subsequently dump the validation NTuple. 
+if __name__=="__main__":
+    from MuonGeoModelTestR4.testGeoModel import setupGeoR4TestCfg, SetupArgParser, executeTest,setupHistSvcCfg
+    parser = SetupArgParser()
+    parser.set_defaults(nEvents = -1)
+
+    args = parser.parse_args()
+    flags, cfg = setupGeoR4TestCfg(args)
+    
+    from xAODMuonSimHitCnv.MuonSimHitCnvCfg import xAODSimHitToMdtMeasCnvAlgCfg
+    from MuonFastDigiTest.MuonFastDigiTesterConfig import MDTFastDigiTesterCfg
+
+    cfg.merge(setupHistSvcCfg(flags,out_file=args.outRootFile,out_stream="MuonFastDigiTest"))
+    cfg.merge(xAODSimHitToMdtMeasCnvAlgCfg(flags))
+    cfg.merge(MDTFastDigiTesterCfg(flags))
+
+    # output spam reduction
+    cfg.getService("AthenaHiveEventLoopMgr").EventPrintoutInterval=500
+
+
+    executeTest(cfg, args.nEvents)
+    

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonFastDigiTest/src/MDTFastDigiTester.cxx

+/*
+  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MDTFastDigiTester.h"
+
+
+namespace MuonValR4 {
+    MDTFastDigiTester::MDTFastDigiTester(const std::string& name,
+                                ISvcLocator* pSvcLocator)
+        : AthHistogramAlgorithm(name, pSvcLocator) {}
+
+    MDTFastDigiTester::~MDTFastDigiTester() {}
+
+    StatusCode MDTFastDigiTester::initialize() {
+        ATH_CHECK(m_inSimHitKey.initialize());
+        ATH_CHECK(m_tree.init(this));
+        ATH_CHECK(m_inDriftCircleKey.initialize());
+        ATH_CHECK(m_idHelperSvc.retrieve());
+        ATH_MSG_DEBUG("Succesfully initialised");
+        return StatusCode::SUCCESS;
+    }
+
+    StatusCode MDTFastDigiTester::finalize() {
+        ATH_CHECK(m_tree.write());
+        return StatusCode::SUCCESS;
+    }
+
+    StatusCode MDTFastDigiTester::execute()  {
+
+        // retrieve the two input collections
+        const EventContext & context = Gaudi::Hive::currentContext();
+        SG::ReadHandle<xAOD::MuonSimHitContainer> readSimHits(m_inSimHitKey,
+                                                              context);
+        SG::ReadHandle<xAOD::MdtDriftCircleContainer> readDriftCircles(
+            m_inDriftCircleKey, context);
+        if (!readSimHits.isValid()){
+            ATH_MSG_FATAL("Failed to read MDT sim hits");
+        }
+        if (!readDriftCircles.isValid()){
+            ATH_MSG_FATAL("Failed to read MDT drift circles");
+        }
+        ATH_MSG_DEBUG("Succesfully retrieved input collections");
+
+        // map the drift circles to identifiers. 
+        // The fast digi should only generate one circle per tube. 
+        std::map<Identifier, const xAOD::MdtDriftCircle*> driftCircleMap{};
+
+        for (const xAOD::MdtDriftCircle* driftCircle : *readDriftCircles) {
+            // need a small type conversion, as the UncalibratedMeasurement::identifier
+            // returns a 32-bit unsigned long, while the Identifier::value_type is 64-bit
+            Identifier driftCircleID{(Identifier::value_type)driftCircle->identifier()};
+            auto empl_res = driftCircleMap.emplace(driftCircleID,driftCircle); 
+            // report any found duplication
+            if (!empl_res.second){
+                ATH_MSG_WARNING("Duplicate drift circle found for identifier "<<m_idHelperSvc->toString(driftCircleID));
+                ATH_MSG_WARNING("  This instance is at r = "<<driftCircle->driftRadius());
+                ATH_MSG_WARNING("  The conflicting instance is at r = "<<empl_res.first->second->driftRadius()); 
+            }
+        }
+        // also track duplicate sim hits for early validation
+        std::map<Identifier,const xAOD::MuonSimHit*> simHitIDsSeen;
+
+        for (const xAOD::MuonSimHit* simHit : *readSimHits) {
+            // The sim hit collection contains non-muon hits, while the fast digi only writes out muon 
+            // deposits. Here remove sim hits that are not expected to be accounted for in the digi. 
+            if (std::abs(simHit->pdgId()) != 13) continue;
+            Identifier ID = simHit->identify();    
+            // check for duplicates 
+            auto empl_res = simHitIDsSeen.emplace(ID,simHit); 
+            // report any found duplication
+            if (!empl_res.second){
+                ATH_MSG_WARNING("Duplicate sim hit found for identifier "<<m_idHelperSvc->toString(ID));
+                ATH_MSG_WARNING("  This instance is at r,z "
+                                    <<simHit->localPosition().perp()
+                                    <<", "<<simHit->localPosition().z()
+                                    <<" from a "<<simHit->pdgId()
+                                    <<" with BC "<<simHit->genParticleLink().barcode()); 
+                ATH_MSG_WARNING("  The conflicting instance is at r,z "
+                                    <<empl_res.first->second->localPosition().perp()
+                                    <<", "<<empl_res.first->second->localPosition().z()
+                                    <<" from a "<<empl_res.first->second->pdgId()
+                                    <<" with BC "<<empl_res.first->second->genParticleLink().barcode()); 
+            }
+            const MdtIdHelper& mdtHelper{m_idHelperSvc->mdtIdHelper()};
+            // populate initial output
+            m_out_stationName = m_idHelperSvc->stationName(ID); 
+            m_out_stationEta = m_idHelperSvc->stationEta(ID); 
+            m_out_stationPhi = m_idHelperSvc->stationPhi(ID); 
+            m_out_multilayer = mdtHelper.multilayer(ID); 
+            m_out_tubeLayer = mdtHelper.tubeLayer(ID); 
+            m_out_tube = mdtHelper.tube(ID); 
+            m_out_simDriftRadius = simHit->localPosition().perp(); 
+            // find the corresponding drift circle
+            auto foundDriftCircle = driftCircleMap.find(ID);
+            if(foundDriftCircle != driftCircleMap.end()){
+                const xAOD::MdtDriftCircle* dt = foundDriftCircle->second; 
+                m_out_digiDriftRadius = dt->driftRadius();
+                m_out_digiDriftRadiusCov = dt->driftRadiusCov();
+                m_out_hasDigi = true;             
+            } 
+            // fill the tree
+            ATH_CHECK(m_tree.fill(context)); 
+        }
+        return StatusCode::SUCCESS;
+    }
+
+}  // namespace MuonValR4

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonFastDigiTest/src/MDTFastDigiTester.h

+/*
+  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MUONFASTDIGITEST_MUONVALR4_MDTFastDigiTester_H
+#define MUONFASTDIGITEST_MUONVALR4_MDTFastDigiTester_H
+
+// Framework includes
+#include "AthenaBaseComps/AthHistogramAlgorithm.h"
+#include "StoreGate/ReadHandleKey.h"
+
+// EDM includes 
+#include "xAODMuonSimHit/MuonSimHitContainer.h"
+#include "xAODMuonPrepData/MdtDriftCircleContainer.h"
+
+// muon includes
+#include "MuonIdHelpers/IMuonIdHelperSvc.h"
+#include "MuonTesterTree/MuonTesterTree.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.
+
+namespace MuonValR4{
+
+  class MDTFastDigiTester : public AthHistogramAlgorithm {
+  public:
+    MDTFastDigiTester(const std::string& name, ISvcLocator* pSvcLocator);
+    virtual ~MDTFastDigiTester() override;
+
+    virtual StatusCode initialize() override;
+    virtual StatusCode execute() override;
+    virtual StatusCode finalize() override;
+
+  private:
+    // MDT sim hits in xAOD format 
+    SG::ReadHandleKey<xAOD::MuonSimHitContainer> m_inSimHitKey {this, "SimHitKey", "xMdtSimHits",
+                                                          "xAOD SimHit collection"};
+    // drift circles in xAOD format 
+    SG::ReadHandleKey<xAOD::MdtDriftCircleContainer> m_inDriftCircleKey{
+    this, "DriftCircleKey", "xAODMdtCircles", "mdt circle container"};
+    
+    ServiceHandle<Muon::IMuonIdHelperSvc> m_idHelperSvc{
+        this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"};
+
+    // output tree - allows to compare the sim and fast-digitised hits
+    MuonVal::MuonTesterTree m_tree{"FastSmearingResults","MuonFastDigiTest"}; 
+    MuonVal::ScalarBranch<int> m_out_stationName{m_tree,"stationName",-1};
+    MuonVal::ScalarBranch<int> m_out_stationEta{m_tree,"stationEta",-1};
+    MuonVal::ScalarBranch<int> m_out_stationPhi{m_tree,"stationPhi",-1};
+    MuonVal::ScalarBranch<int> m_out_multilayer{m_tree,"multilayer",-1};
+    MuonVal::ScalarBranch<int> m_out_tubeLayer{m_tree,"tubeLayer",-1};
+    MuonVal::ScalarBranch<int> m_out_tube{m_tree,"tube",-1};    
+    MuonVal::ScalarBranch<double> m_out_simDriftRadius{m_tree, "sim_driftRadius", 0.0}; 
+    MuonVal::ScalarBranch<bool> m_out_hasDigi{m_tree, "hasDigi", false}; 
+    MuonVal::ScalarBranch<double> m_out_digiDriftRadius{m_tree, "digi_driftRadius", 0.0}; 
+    MuonVal::ScalarBranch<double> m_out_digiDriftRadiusCov{m_tree, "digi_driftRadiusCov", 0.0}; 
+  };
+}
+
+#endif // MUONFASTDIGITEST_MUONVALR4_MDTFastDigiTester_H

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonFastDigiTest/src/components/MuonFastDigiTest_entries.cxx

+/*
+  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "../MDTFastDigiTester.h"
+
+DECLARE_COMPONENT(MuonValR4::MDTFastDigiTester)