ATLAS Offline Software
MuonSDOAnalysis.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "MuonSDOAnalysis.h"
6 #include "StoreGate/ReadHandle.h"
7 #include "GaudiKernel/SystemOfUnits.h"
8 
9 #include <algorithm>
10 #include <format>
11 
12 #include <TH1F.h>
13 #include <TH2F.h>
14 namespace MuonVal{
15  using namespace Muon::MuonStationIndex;
16 
19  const std::string& basePath,
20  const TechIdx_t techIdx,
21  const int stIdent):
22  identifier{stIdent}{
23 
24  auto registerHisto = [&](TH1* h) {
25  h->Sumw2();
26  h->StatOverflows();
27  std::string stName = "Inclusive";
28  if (stIdent >= 0) {
29  switch (techIdx) {
30  case TechIdx_t::MDT:
31  case TechIdx_t::RPC: {
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 
53  h_sdoWord = registerHisto(new TH1F("sdoWord", "sdoWord", 100, 0, 10));
54  h_eventIndex = registerHisto(new TH1F("eventIndex", "Event index (SDO)", 100, 0, 1000));
55 
56  switch (techIdx) {
57  case TechIdx_t::TGC: {
58  h_stationPhi = registerHisto(new TH1F("stationPhi", "StationPhi;stationPhi",48, 0.5, 48.5));
59  h_stationEta = registerHisto(new TH1F("stationEta", "StationEta;stationEta",11, -5.5, 5.5));
60  constexpr double r = 12. * Gaudi::Units::m;
61  constexpr double z = 14 * Gaudi::Units::m;
62  h_xyPos = registerHisto(new TH2F("hitXY", "hitXY;x [mm]; y[mm];entries", 240, -r,r, 240,-r,r));
63  h_rzPos = registerHisto(new TH2F("hitRZ", "hitRz;r [mm]; z[mm];entries", 1200, 0, r, 2200,-z,z));
64  break;
65  } case TechIdx_t::RPC: {
66  h_stationPhi = registerHisto(new TH1F("stationPhi", "StationPhi;stationPhi",8, 0.5, 8.5));
67  h_stationEta = registerHisto(new TH1F("stationEta", "StationEta;stationEta",17, -8.5, 8.5));
68  constexpr double r = 12. * Gaudi::Units::m;
69  constexpr double z = 12 * Gaudi::Units::m;
70  h_xyPos = registerHisto(new TH2F("hitXY", "hitXY;x [mm]; y[mm];entries", 240, -r,r, 240,-r,r));
71  h_rzPos = registerHisto(new TH2F("hitRZ", "hitRz;r [mm]; z[mm];entries", 1200, 0, r, 600,-z,z));
72 
73  break;
74  } case TechIdx_t::MDT: {
75  h_stationPhi = registerHisto(new TH1F("stationPhi", "StationPhi;stationPhi",8, 0.5, 8.5));
76  h_stationEta = registerHisto(new TH1F("stationEta", "StationEta;stationEta",17, -8.5, 8.5));
77  constexpr double r = 12. * Gaudi::Units::m;
78  constexpr double z = 22 * Gaudi::Units::m;
79  h_xyPos = registerHisto(new TH2F("hitXY", "hitXY;x [mm]; y[mm];entries", 2400, -r,r, 2400,-r,r));
80  h_rzPos = registerHisto(new TH2F("hitRZ", "hitRz;r [mm]; z[mm];entries", 1200, 0, r, 2200,-z,z));
81 
82  break;
83  } case TechIdx_t::MM: {
84  h_stationPhi = registerHisto(new TH1F("stationPhi", "StationPhi;stationPhi",8, 0.5, 8.5));
85  h_stationEta = registerHisto(new TH1F("stationEta", "StationEta;stationEta",5, -2.5, 2.5));
86  constexpr double r = 5 * Gaudi::Units::m;
87  constexpr double z = 8. *Gaudi::Units::m;
88  h_xyPos = registerHisto(new TH2F("hitXY", "hitXY;x [mm]; y[mm];entries", 1000, -r,r, 1000,-r,r));
89  h_rzPos = registerHisto(new TH2F("hitRZ", "hitRz;r [mm]; z[mm];entries", 500, 0, r, 160, -z, z));
90 
91  break;
92  } case TechIdx_t::STGC: {
93  h_stationPhi = registerHisto(new TH1F("stationPhi", "StationPhi;stationPhi",8, 0.5, 8.5));
94  h_stationEta = registerHisto(new TH1F("stationEta", "StationEta;stationEta",7, -3.5, 3.5));
95  constexpr double r = 5 * Gaudi::Units::m;
96  constexpr double z = 8. *Gaudi::Units::m;
97  h_xyPos = registerHisto(new TH2F("hitXY", "hitXY;x [mm]; y[mm];entries", 1000, -r,r, 1000,-r,r));
98  h_rzPos = registerHisto(new TH2F("hitRZ", "hitRz;r [mm]; z[mm];entries", 500, 0, r, 160, -z, z));
99  break;
100  } default:
101  break;
102  }
103 
104  h_radius = registerHisto(new TH1F("h_radius", "Radius (SDO)", 100, -15, 15));
105  h_localZ = registerHisto(new TH1F("h_localZ", "Local z-position (SDO)", 100, -3200, 3200));
106  }
107 
109 
111  ATH_CHECK(m_idHelperSvc.retrieve());
112  std::set<int> histIds{-1};
113  TechIdx_t techIdx{static_cast<TechIdx_t>(1*m_techIdx)};
114  using namespace Muon::MuonStationIndex;
115  switch (techIdx) {
116  case TechIdx_t::MDT: {
117  for (int s = 0 ; s < toInt(StIdx_t::StIndexMax); ++s){
118  histIds.insert(s);
119  }
120  break;
121  } case TechIdx_t::RPC: {
122  for (int s = 0 ; s < toInt(StIdx_t::BE); ++s){
123  histIds.insert(s);
124  }
125  break;
126  } case TechIdx_t::TGC: {
127  const auto& idHelper{m_idHelperSvc->tgcIdHelper()};
128  std::for_each(idHelper.module_begin(), idHelper.module_end(),
129  [&histIds, &idHelper](const Identifier& id) {
130  histIds.insert(idHelper.stationName(id));
131  });
132  break;
133  } case TechIdx_t::STGC: {
134  const auto& idHelper{m_idHelperSvc->stgcIdHelper()};
135  std::for_each(idHelper.module_begin(), idHelper.module_end(),
136  [&histIds, &idHelper](const Identifier& id) {
137  histIds.insert(idHelper.stationName(id));
138  });
139  break;
140  } case TechIdx_t::MM: {
141  const auto& idHelper{m_idHelperSvc->mmIdHelper()};
142  std::for_each(idHelper.module_begin(), idHelper.module_end(),
143  [&histIds, &idHelper](const Identifier& id) {
144  histIds.insert(idHelper.stationName(id));
145  });
146  break;
147  } default:
148  ATH_MSG_FATAL("Invalid technology index parsed "<<m_techIdx);
149  return StatusCode::FAILURE;
150  }
151  for (int id : histIds) {
152  m_histos.emplace_back(m_idHelperSvc.get(), histSvc(), m_path, techIdx, id);
153  }
154  return StatusCode::SUCCESS;
155  }
156 
158  const EventContext& ctx{Gaudi::Hive::currentContext()};
159  const MuonSimDataCollection* sdoCont{nullptr};
160  ATH_CHECK(SG::get(sdoCont, m_inputKey, ctx));
161 
162  for (const auto&[sdoID, sdo]: *sdoCont) {
163  std::vector<unsigned> fillMe{};
164  const TechIdx_t techIdx = m_idHelperSvc->technologyIndex(sdoID);
165 
166  const int stIdx = (techIdx == TechIdx_t::TGC || techIdx == TechIdx_t::STGC ||
167  techIdx == TechIdx_t::MM ? m_idHelperSvc->stationName(sdoID)
168  : toInt(m_idHelperSvc->stationIndex(sdoID)));
169  const int stationEta = m_idHelperSvc->stationEta(sdoID);
170  const int stationPhi = m_idHelperSvc->stationPhi(sdoID);
171 
172 
173  for (unsigned h = 0 ; h < m_histos.size(); ++h) {
174  if (m_histos[h].identifier == -1 ||
175  m_histos[h].identifier == stIdx) {
176  fillMe.push_back(h);
177  }
178  }
179 
180  const int sdoWord = sdo.word();
181  const Amg::Vector3D gPos = sdo.globalPosition();
182  const float xPos = gPos.x();
183  const float yPos = gPos.y();
184  const float zPos = gPos.z();
185  for (unsigned h : fillMe) {
186  m_histos[h].h_sdoWord->Fill(sdoWord);
187  m_histos[h].h_xyPos->Fill(xPos, yPos);
188  m_histos[h].h_rzPos->Fill(std::hypot(xPos, yPos), zPos);
189  m_histos[h].h_stationPhi->Fill(stationPhi);
190  m_histos[h].h_stationEta->Fill(stationEta);
191 
192  }
193  for (const auto& [particleLink, data]: sdo.getdeposits()) {
194  const int eventIx = particleLink.eventIndex();
195  const float radius = data.firstEntry();
196  const float localZ = data.secondEntry();
197  for (unsigned h : fillMe) {
198  m_histos[h].h_eventIndex->Fill(eventIx);
199  m_histos[h].h_radius->Fill(radius);
200  m_histos[h].h_localZ->Fill(localZ);
201  }
202  }
203  }
204  return StatusCode::SUCCESS;
205  }
206 }
207 
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
beamspotman.r
def r
Definition: beamspotman.py:672
data
char data[hepevt_bytes_allocation_ATLAS]
Definition: HepEvt.cxx:11
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
vtune_athena.format
format
Definition: vtune_athena.py:14
Muon::MuonStationIndex
Definition: MuonStationIndex.h:13
Muon::MuonStationIndex::TechnologyIndex::RPC
@ RPC
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
MuonSDOAnalysis.h
MuonVal::MuonSDOAnalysis::m_inputKey
SG::ReadHandleKey< MuonSimDataCollection > m_inputKey
Input read handle key.
Definition: MuonSDOAnalysis.h:27
xAOD::identifier
identifier
Definition: UncalibratedMeasurement_v1.cxx:15
MuonVal::MuonSDOAnalysis::initialize
virtual StatusCode initialize() override final
Definition: MuonSDOAnalysis.cxx:108
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
MuonVal::MuonSDOAnalysis::execute
virtual StatusCode execute() override final
Definition: MuonSDOAnalysis.cxx:157
MuonVal::MuonSDOAnalysis::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Service handle of the IdHelperSvc.
Definition: MuonSDOAnalysis.h:29
Muon::MuonStationIndex::TechnologyIndex::MM
@ MM
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
extractSporadic.h
list h
Definition: extractSporadic.py:96
MuonSimDataCollection
Definition: MuonSimDataCollection.h:21
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
Muon::MuonStationIndex::StIndex
StIndex
enum to classify the different station layers in the muon spectrometer
Definition: MuonStationIndex.h:23
MuonVal
Class to store array like branches into the n-tuples.
Definition: HitValAlg.cxx:19
Muon::MuonStationIndex::StIndex::BE
@ BE
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
MuonVal::MuonSDOAnalysis::m_techIdx
Gaudi::Property< int > m_techIdx
Definition: MuonSDOAnalysis.h:36
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
LVL1::gFEXPos
Definition: gFexPos.h:11
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
MuonVal::MuonSDOAnalysis::m_path
Gaudi::Property< std::string > m_path
Definition: MuonSDOAnalysis.h:31
MuonVal::MuonSDOAnalysis::HistoSet::HistoSet
HistoSet()=default
Default constructor.
ReadHandle.h
Handle class for reading from StoreGate.
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
Muon::MuonStationIndex::toInt
constexpr int toInt(const EnumType enumVal)
Definition: MuonStationIndex.h:61
ServiceHandle< ITHistSvc >
MuonVal::MuonSDOAnalysis::m_histos
std::vector< HistoSet > m_histos
Definition: MuonSDOAnalysis.h:79
Identifier
Definition: IdentifierFieldParser.cxx:14