ATLAS Offline Software
xMuonHitAnalysis.cxx
Go to the documentation of this file.
1 
2 /*
3  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
4 */
5 #include "xMuonHitAnalysis.h"
6 
7 #include "StoreGate/ReadHandle.h"
8 #include "xAODMuonViews/ChamberViewer.h"
9 #include "GaudiKernel/SystemOfUnits.h"
10 #include "GeoPrimitives/GeoPrimitivesHelpers.h"
11 
12 #include <format>
13 #include <TH1F.h>
14 #include <TH2F.h>
15 namespace MuonValR4{
16  xMuonHitAnalysis::HistoSet::HistoSet(const Muon::IMuonIdHelperSvc* idHelperSvc,
17  const ServiceHandle<ITHistSvc>& histSvc,
18  const std::string& basePath,
19  const TechIdx_t techIdx,
20  const int stIdent):
21  identifier{stIdent} {
22 
23  auto registerHisto = [&](TH1* h) {
24  h->Sumw2();
25  h->StatOverflows();
26 
27  std::string stName = "Inclusive";
28  if (stIdent >= 0) {
29  switch (techIdx) {
30  case TechIdx_t::MDT:
31  case TechIdx_t::RPC: {
32  stName = Muon::MuonStationIndex::stName(static_cast<StIdx_t>(identifier));
33  break;
34  } case TechIdx_t::TGC: {
35  stName = idHelperSvc->tgcIdHelper().stationNameString(stIdent);
36  break;
37  } case TechIdx_t::STGC: {
38  stName = idHelperSvc->stgcIdHelper().stationNameString(stIdent);
39  break;
40  } case TechIdx_t::MM: {
41  stName = idHelperSvc->mmIdHelper().stationNameString(stIdent);
42  break;
43  } default:
44  break;
45  }
46  }
47  const std::string histPath = std::format("/{}/Station_{}/{}",
48  basePath, stName, h->GetName());
49  histSvc->regHist(histPath, h).ignore();
50  return h;
51  };
52  h_dirTheta = registerHisto(new TH1F("dirTheta", "local direction; #theta [deg]", 180, 0, 180));
53  h_dirPhi = registerHisto(new TH1F("dirPhi", "local direction; #phi [deg]", 360, -180, 180));
54  h_energy = registerHisto(new TH1F("energy", "kinetic energy; e[GeV]", 150, 0, 150));
55  h_deposit = registerHisto(new TH1F("energyDeposit", "energy deposit; e[GeV]", 20, 0, 2));
56  h_pdgId = registerHisto(new TH1F("pdgId", "Particle data ID; pdgId",31, -15.5, 15.5));
57 
58  switch (techIdx) {
59  case TechIdx_t::TGC: {
60  h_stationPhi = registerHisto(new TH1F("stationPhi", "StationPhi;stationPhi",48, 0.5, 48.5));
61  h_stationEta = registerHisto(new TH1F("stationEta", "StationEta;stationEta",11, -5.5, 5.5));
62  h_localHitPos = registerHisto(new TH2F("localHitPosXY", "hit position;x [mm]; y[mm]",
63  240, -1200, 1200, 240, -1200, 1200));
64  h_localHitZ = registerHisto(new TH1F("localHitPosZ", "position in gas gap;z [mm]", 20, -5, 5));
65 
66  constexpr double r = 12. * Gaudi::Units::m;
67  constexpr double z = 14 * Gaudi::Units::m;
68  h_globHitXY = registerHisto(new TH2F("hitXY", "hitXY;x [mm]; y[mm];entries", 240, -r,r, 240,-r,r));
69  h_globHitRZ = registerHisto(new TH2F("hitRZ", "hitRz;r [mm]; z[mm];entries", 1200, 0, r, 2200,-z,z));
70  break;
71  } case TechIdx_t::RPC: {
72  h_stationPhi = registerHisto(new TH1F("stationPhi", "StationPhi;stationPhi",8, 0.5, 8.5));
73  h_stationEta = registerHisto(new TH1F("stationEta", "StationEta;stationEta",17, -8.5, 8.5));
74  h_localHitPos = registerHisto(new TH2F("localHitPosXY", "hit position;x [mm]; y[mm]",
75  70, -1400, 1400, 70, -1400, 1400));
76  h_localHitZ = registerHisto(new TH1F("localHitPosZ", "position in gas gap;z [mm]", 11, -0.5, 0.5));
77  constexpr double r = 12. * Gaudi::Units::m;
78  constexpr double z = 12 * Gaudi::Units::m;
79  h_globHitXY = registerHisto(new TH2F("hitXY", "hitXY;x [mm]; y[mm];entries", 240, -r,r, 240,-r,r));
80  h_globHitRZ = registerHisto(new TH2F("hitRZ", "hitRz;r [mm]; z[mm];entries", 1200, 0, r, 600,-z,z));
81 
82  break;
83  } case TechIdx_t::MDT: {
84  h_stationPhi = registerHisto(new TH1F("stationPhi", "StationPhi;stationPhi",8, 0.5, 8.5));
85  h_stationEta = registerHisto(new TH1F("stationEta", "StationEta;stationEta",17, -8.5, 8.5));
86  h_localHitPos = registerHisto(new TH1F("localHitPosR", "hit position;r [mm]", 61, -15, 15));
87  h_localHitZ = registerHisto(new TH1F("localHitPosZ", "position in gas gap;z [mm]", 310, -3100, 3100));
88  constexpr double r = 12. * Gaudi::Units::m;
89  constexpr double z = 22 * Gaudi::Units::m;
90  h_globHitXY = registerHisto(new TH2F("hitXY", "hitXY;x [mm]; y[mm];entries", 2400, -r,r, 2400,-r,r));
91  h_globHitRZ = registerHisto(new TH2F("hitRZ", "hitRz;r [mm]; z[mm];entries", 1200, 0, r, 2200,-z,z));
92 
93  break;
94  } case TechIdx_t::MM: {
95  h_stationPhi = registerHisto(new TH1F("stationPhi", "StationPhi;stationPhi",8, 0.5, 8.5));
96  h_stationEta = registerHisto(new TH1F("stationEta", "StationEta;stationEta",5, -2.5, 2.5));
97  h_localHitPos = registerHisto(new TH2F("localHitPosXY", "hit position;x [mm]; y[mm]",
98  480, -1200, 1200, 480, -2400, 2400));
99  h_localHitZ = registerHisto(new TH1F("localHitPosZ", "position in gas gap;z [mm]", 11, -0.5, 0.5));
100  constexpr double r = 5 * Gaudi::Units::m;
101  constexpr double z = 8. *Gaudi::Units::m;
102  h_globHitXY = registerHisto(new TH2F("hitXY", "hitXY;x [mm]; y[mm];entries", 1000, -r,r, 1000,-r,r));
103  h_globHitRZ = registerHisto(new TH2F("hitRZ", "hitRz;r [mm]; z[mm];entries", 500, 0, r, 160, -z, z));
104 
105  break;
106  } case TechIdx_t::STGC: {
107  h_stationPhi = registerHisto(new TH1F("stationPhi", "StationPhi;stationPhi",8, 0.5, 8.5));
108  h_stationEta = registerHisto(new TH1F("stationEta", "StationEta;stationEta",7, -3.5, 3.5));
109  h_localHitPos = registerHisto(new TH2F("localHitPosXY", "hit position;x [mm]; y[mm]",
110  440, -1100, 1100, 440, -1100, 1100));
111  h_localHitZ = registerHisto(new TH1F("localHitPosZ", "position in gas gap;z [mm]", 11, -0.5, 0.5));
112  constexpr double r = 5 * Gaudi::Units::m;
113  constexpr double z = 8. *Gaudi::Units::m;
114  h_globHitXY = registerHisto(new TH2F("hitXY", "hitXY;x [mm]; y[mm];entries", 1000, -r,r, 1000,-r,r));
115  h_globHitRZ = registerHisto(new TH2F("hitRZ", "hitRz;r [mm]; z[mm];entries", 500, 0, r, 160, -z, z));
116  break;
117  } default:
118  break;
119  }
120  }
121 
122 
123  StatusCode xMuonHitAnalysis::initialize(){
124  ATH_CHECK(m_inputKey.initialize());
125  ATH_CHECK(m_geoCtxKey.initialize());
126  ATH_CHECK(m_idHelperSvc.retrieve());
127  ATH_CHECK(detStore()->retrieve(m_detMgr));
128  std::set<int> histIds{-1};
129  TechIdx_t techIdx{static_cast<TechIdx_t>(1*m_techIdx)};
130  using namespace Muon::MuonStationIndex;
131  switch (techIdx) {
132  case TechIdx_t::MDT: {
133  for (int s = 0 ; s < toInt(StIdx_t::StIndexMax); ++s){
134  histIds.insert(s);
135  }
136  break;
137  } case TechIdx_t::RPC: {
138  for (int s = 0 ; s < toInt(StIdx_t::BE); ++s){
139  histIds.insert(s);
140  }
141  break;
142  } case TechIdx_t::TGC: {
143  const auto& idHelper{m_idHelperSvc->tgcIdHelper()};
144  std::for_each(idHelper.module_begin(), idHelper.module_end(),
145  [&histIds, &idHelper](const Identifier& id) {
146  histIds.insert(idHelper.stationName(id));
147  });
148  break;
149  } case TechIdx_t::STGC: {
150  const auto& idHelper{m_idHelperSvc->stgcIdHelper()};
151  std::for_each(idHelper.module_begin(), idHelper.module_end(),
152  [&histIds, &idHelper](const Identifier& id) {
153  histIds.insert(idHelper.stationName(id));
154  });
155  break;
156  } case TechIdx_t::MM: {
157  const auto& idHelper{m_idHelperSvc->mmIdHelper()};
158  std::for_each(idHelper.module_begin(), idHelper.module_end(),
159  [&histIds, &idHelper](const Identifier& id) {
160  histIds.insert(idHelper.stationName(id));
161  });
162  break;
163  } default:
164  ATH_MSG_FATAL("Invalid technology index parsed "<<m_techIdx);
165  return StatusCode::FAILURE;
166  }
167  for (int id : histIds) {
168  m_histos.emplace_back(m_idHelperSvc.get(), histSvc(), m_path, techIdx, id);
169  }
170  return StatusCode::SUCCESS;
171  }
172  StatusCode xMuonHitAnalysis::execute() {
173  const EventContext& ctx{Gaudi::Hive::currentContext()};
174  const xAOD::MuonSimHitContainer* simHits{nullptr};
175  const ActsGeometryContext* gctx{nullptr};
176  ATH_CHECK(SG::get(simHits, m_inputKey, ctx));
177  if (simHits->empty()) {
178  ATH_MSG_DEBUG("No hits recorded in the event "<<m_inputKey.fullKey());
179  return StatusCode::SUCCESS;
180  }
181  ATH_CHECK(SG::get(gctx, m_geoCtxKey, ctx));
182 
183  xAOD::ChamberViewer viewer{*simHits, m_idHelperSvc.get(),
184  xAOD::ChamberView::Mode::Chamber};
185  do{
186  const Identifier refId{viewer.at(0)->identify()};
187  const TechIdx_t techIdx = m_idHelperSvc->technologyIndex(refId);
188 
189  const int stIdx = (techIdx == TechIdx_t::TGC || techIdx == TechIdx_t::STGC ||
190  techIdx == TechIdx_t::MM ? m_idHelperSvc->stationName(refId)
191  : toInt(m_idHelperSvc->stationIndex(refId)));
192  const int stationEta = m_idHelperSvc->stationEta(refId);
193  const int stationPhi = m_idHelperSvc->stationPhi(refId);
194 
195  std::vector<unsigned> fillMe{};
196  for (unsigned h = 0; h<m_histos.size(); ++h) {
197  if (m_histos[h].identifier == -1 ||
198  m_histos[h].identifier == stIdx){
199  fillMe.push_back(h);
200  }
201  }
202  for (const xAOD::MuonSimHit* hit : viewer) {
203  const Identifier& hitId{hit->identify()};
204  const MuonGMR4::MuonReadoutElement* re = m_detMgr->getReadoutElement(hitId);
205  const double theta = hit->localDirection().theta();
206  const double phi = hit->localDirection().phi();
207  double signedR{0.};
208  IdentifierHash layHash{re->layerHash(hitId)};
209  if (techIdx ==TechIdx_t::MDT) {
210  signedR = Amg::lineDistance<3>(Amg::Vector3D::Zero(),
211  Amg::Vector3D::UnitZ(),
212  xAOD::toEigen(hit->localPosition()),
213  xAOD::toEigen(hit->localDirection()));
214  layHash = re->measurementHash(hitId);
215  }
216  const Amg::Vector3D globPos{re->localToGlobalTrans(*gctx,layHash)*
217  xAOD::toEigen(hit->localPosition())};
218 
219  for (unsigned h : fillMe) {
220  m_histos[h].h_dirTheta->Fill(theta / Gaudi::Units::deg);
221  m_histos[h].h_dirPhi->Fill(phi / Gaudi::Units::deg);
222  m_histos[h].h_energy->Fill(hit->kineticEnergy()/ Gaudi::Units::GeV);
223  m_histos[h].h_deposit->Fill(hit->energyDeposit() / Gaudi::Units::GeV);
224  m_histos[h].h_pdgId->Fill(hit->pdgId());
225  m_histos[h].h_stationPhi->Fill(stationPhi);
226  m_histos[h].h_stationEta->Fill(stationEta);
227  if (techIdx == TechIdx_t::MDT) {
228  m_histos[h].h_localHitPos->Fill(signedR);
229  } else {
230  m_histos[h].h_localHitPos->Fill(hit->localPosition().x(),
231  hit->localPosition().y());
232  }
233  m_histos[h].h_localHitZ->Fill(hit->localPosition().z());
234  m_histos[h].h_globHitXY->Fill(globPos.x(), globPos.y());
235  m_histos[h].h_globHitRZ->Fill(globPos.perp(), globPos.y());
236  }
237  }
238  } while(viewer.next());
239  return StatusCode::SUCCESS;
240  }
241 
242 }
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
AthHistogramAlgorithm::histSvc
const ServiceHandle< ITHistSvc > & histSvc() const
The standard THistSvc (for writing histograms and TTrees and more to a root file) Returns (kind of) a...
Definition: AthHistogramAlgorithm.h:113
MuonValR4::xMuonHitAnalysis::initialize
virtual StatusCode initialize() override final
Definition: xMuonHitAnalysis.cxx:123
beamspotman.r
def r
Definition: beamspotman.py:672
MuonValR4::xMuonHitAnalysis::m_inputKey
SG::ReadHandleKey< xAOD::MuonSimHitContainer > m_inputKey
Definition: xMuonHitAnalysis.h:24
xAOD::MuonSimHit_v1
Definition: MuonSimHit_v1.h:18
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:18
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
vtune_athena.format
format
Definition: vtune_athena.py:14
Muon::MuonStationIndex
Definition: MuonStationIndex.h:13
Muon::MuonStationIndex::TechnologyIndex::RPC
@ RPC
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:75
Muon::MuonStationIndex::stName
const std::string & stName(StIndex index)
convert StIndex into a string
Definition: MuonStationIndex.cxx:104
Muon::MuonStationIndex::TechnologyIndex
TechnologyIndex
enum to classify the different layers in the muon spectrometer
Definition: MuonStationIndex.h:54
MuonGMR4::MuonReadoutElement
The MuonReadoutElement is an abstract class representing the geometry representing the muon detector.
Definition: MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/MuonReadoutElement.h:38
deg
#define deg
Definition: SbPolyhedron.cxx:17
xAOD::identifier
identifier
Definition: UncalibratedMeasurement_v1.cxx:15
python.TrigEgammaMonitorHelper.TH2F
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:45
Muon::MuonStationIndex::TechnologyIndex::MDT
@ MDT
AthCommonDataStore< AthCommonMsg< Algorithm > >::detStore
const ServiceHandle< StoreGateSvc > & detStore() const
The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:95
MuonValR4::xMuonHitAnalysis::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Service handle of the IdHelperSvc.
Definition: xMuonHitAnalysis.h:28
xAOD::ChamberView::Mode::Chamber
@ Chamber
View ends if the moduleHash changes.
xAOD::ChamberViewer
Definition: ChamberViewer.h:59
MuonValR4::xMuonHitAnalysis::HistoSet::HistoSet
HistoSet(const Muon::IMuonIdHelperSvc *idHelperSvc, const ServiceHandle< ITHistSvc > &histSvc, const std::string &basePath, const TechIdx_t techIdx, const int stIdent)
Helper struct to define histograms spectrometer region.
Definition: xMuonHitAnalysis.cxx:16
Muon::MuonStationIndex::TechnologyIndex::MM
@ MM
MuonValR4::xMuonHitAnalysis::m_geoCtxKey
SG::ReadHandleKey< ActsGeometryContext > m_geoCtxKey
Definition: xMuonHitAnalysis.h:26
Muon::MuonStationIndex::TechnologyIndex::TGC
@ TGC
z
#define z
SG::get
const T * get(const ReadCondHandleKey< T > &key, const EventContext &ctx)
Convenience function to retrieve an object given a ReadCondHandleKey.
Definition: ReadCondHandle.h:287
h
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
extractSporadic.h
list h
Definition: extractSporadic.py:96
MuonSegmentReaderConfig.histSvc
histSvc
Definition: MuonSegmentReaderConfig.py:96
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
ActsGeometryContext
Include the GeoPrimitives which need to be put first.
Definition: ActsGeometryContext.h:27
MuonValR4::xMuonHitAnalysis::execute
virtual StatusCode execute() override final
Definition: xMuonHitAnalysis.cxx:172
MuonValR4
Lightweight algorithm to read xAOD MDT sim hits and (fast-digitised) drift circles from SG and fill a...
Definition: IPatternVisualizationTool.h:23
Muon::MuonStationIndex::StIndex
StIndex
enum to classify the different station layers in the muon spectrometer
Definition: MuonStationIndex.h:23
MuonValR4::xMuonHitAnalysis::m_path
Gaudi::Property< std::string > m_path
Definition: xMuonHitAnalysis.h:34
Muon::MuonStationIndex::StIndex::BE
@ BE
MuonValR4::xMuonHitAnalysis::m_techIdx
Gaudi::Property< int > m_techIdx
Definition: xMuonHitAnalysis.h:39
ChamberViewer.h
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
xMuonHitAnalysis.h
h
GeoPrimitivesHelpers.h
re
const boost::regex re(r_e)
Muon::MuonStationIndex::TechnologyIndex::STGC
@ STGC
Muon::IMuonIdHelperSvc
Interface for Helper service that creates muon Identifiers and can be used to print Identifiers.
Definition: IMuonIdHelperSvc.h:27
python.SystemOfUnits.s
float s
Definition: SystemOfUnits.py:147
MuonValR4::xMuonHitAnalysis::m_detMgr
const MuonGMR4::MuonDetectorManager * m_detMgr
MuonSimHit.
Definition: xMuonHitAnalysis.h:32
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
ReadHandle.h
Handle class for reading from StoreGate.
MuonValR4::xMuonHitAnalysis::m_histos
std::vector< HistoSet > m_histos
Definition: xMuonHitAnalysis.h:79
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
Muon::MuonStationIndex::toInt
constexpr int toInt(const EnumType enumVal)
Definition: MuonStationIndex.h:61
MuonGMR4::MuonDetectorManager::getReadoutElement
const MuonReadoutElement * getReadoutElement(const Identifier &id) const
Returns a generic Muon readout element.
ServiceHandle< ITHistSvc >
Identifier
Definition: IdentifierFieldParser.cxx:14
Generated on Mon Sep 29 2025 21:24:03 for ATLAS Offline Software by   doxygen 1.8.18