20 #include "GaudiKernel/PhysicalConstants.h"
22 #include <unordered_map>
29 using namespace SegmentFit;
34 if (m_readKeys.empty()){
35 ATH_MSG_ERROR(
"No simulated hit containers have been parsed to build the segments from ");
36 return StatusCode::FAILURE;
38 ATH_CHECK(m_mdtCalibKey.initialize(m_idHelperSvc->hasMDT()));
39 ATH_CHECK(m_nswUncertKey.initialize(m_idHelperSvc->hasMM() || m_idHelperSvc->hasSTGC()));
50 return StatusCode::SUCCESS;
54 switch (
const auto techIdx = m_idHelperSvc->technologyIndex(hitId)) {
58 if (!
SG::get(calibCont, m_mdtCalibKey, ctx).isSuccess()) {
61 const auto& rtCalib{calibCont->getCalibData(hitId, msgStream())->rtRelation};
62 const double driftTime = rtCalib->tr()->driftTime(hit.
localPosition().perp()).value_or(rtCalib->tr()->maxRadius());
63 return rtCalib->rtRes()->resolution(
driftTime);
65 const auto*
re = m_detMgr->getRpcReadoutElement(hitId);
66 return re->stripEtaPitch() / std::sqrt(12.);
68 const auto*
re = m_detMgr->getTgcReadoutElement(hitId);
69 const IdentifierHash measHash =
re->measurementHash(hitId);
70 const auto& design =
re->wireGangLayout(measHash);
71 return design.stripPitch() / std::sqrt(12.);
75 if (!
SG::get(errorCalibDB, m_nswUncertKey, ctx).isSuccess()) {
81 if (techIdx ==
STGC) {
82 errorCalibInput.clusterAuthor = 3;
84 errorCalibInput.clusterAuthor=66;
86 return errorCalibDB->clusterUncertainty(errorCalibInput);
96 return link->momentum().perp();
116 ATH_MSG_VERBOSE(
"Assemble segments from "<<simHits.size()<<
" background hits.");
117 std::vector<char> alreadyUsed(simHits.size(), 0);
118 for (std::size_t
h = 0 ;
h < simHits.size(); ++
h) {
119 if (alreadyUsed[
h]) {
122 const auto& [refHit, refPos, refDir] = simHits[
h];
124 ATH_MSG_VERBOSE(
"Try to find other hits on trajectory given by "<<m_idHelperSvc->toString(refHit->identify())
126 <<
", pt: "<<muonPt(*refHit, locToGlob.linear()* refDir)/
Gaudi::Units::GeV <<
" [GeV].");
127 std::vector<std::size_t> indicesOnSeg{
h};
128 for (std::size_t
h1 =
h +1;
h1 < simHits.size(); ++
h1) {
129 if (alreadyUsed[
h1]) {
132 const auto&[testHit, testPos, testDir] = simHits[
h1];
135 if (std::abs(refHit->kineticEnergy() - testHit->kineticEnergy()) > m_pileUpHitELoss) {
141 const double angleDir =
Amg::angle(testDir, refDir);
143 ATH_MSG_VERBOSE(
"Test hit "<<m_idHelperSvc->toString(testHit->identify())<<
", "
148 if (std::abs(angleDir) > m_pileUpHitAngleCone || hitSep > m_pileUpHitDistance) {
151 indicesOnSeg.push_back(
h1);
157 globPos, globDir).value_or(0.) * globDir;
158 ATH_MSG_VERBOSE(
"Found "<<indicesOnSeg.size()<<
" matching hits. Closest perigee "
159 <<
Amg::toString(perigee)<<
", "<<m_idCylinderR<<
", "<<m_idCylinderHalfZ);
162 std::back_inserter(hitsOnSeg),
163 [&simHits](
const auto idx) {
166 if (perigee.perp() > m_idCylinderR || std::abs(perigee.z()) > m_idCylinderHalfZ ||
167 constructSegmentFromHits(ctx, locToGlob, hitsOnSeg,
out)){
168 std::ranges::for_each(indicesOnSeg,
169 [&alreadyUsed](
const auto idx) {
170 alreadyUsed[
idx] = 1;
182 auto refHit_itr = std::ranges::min_element(simHits,
184 return std::abs(std::get<1>(
a).
z()) <std::abs(std::get<1>(
b).
z());
187 const auto [simHit, localPos, chamberDir] = (*refHit_itr);
188 ATH_MSG_VERBOSE(
"Create segement from hit: "<<m_idHelperSvc->toString(simHit->identify())<<
189 " pdgId: "<<simHit->pdgId()<<
", energy: "<<simHit->kineticEnergy()
190 <<
", genParticle: "<<simHit->genParticleLink().cptr());
194 const double distance = Amg::intersect<3>(localPos, chamberDir, Amg::Vector3D::UnitZ(), 0.).value_or(0.);
198 const Amg::Vector3D globDir = locToGlob.linear() * chamberDir;
200 unsigned nMdt{0}, nRpcEta{0}, nRpcPhi{0}, nTgcEta{0}, nTgcPhi{0};
201 unsigned nMm{0}, nStgcEta{0}, nStgcPhi{0};
204 for (
const auto& [assocMe,
pos,
dir] : simHits) {
213 if (castRE->nEtaStrips()) ++nRpcEta;
214 if (castRE->nPhiStrips()) ++nRpcPhi;
218 const IdentifierHash gapHash = assocRE->
measurementHash(assocMe->identify());
219 if (castRE->numStrips(gapHash)){
222 if (castRE->numWireGangs(gapHash)) {
234 ATH_MSG_WARNING(
"Csc are not defined "<<m_idHelperSvc->toString(simHit->identify()));
236 ATH_MSG_VERBOSE(
"Associate hit "<<m_idHelperSvc->toString(assocMe->identify())
237 <<
" pdgId: "<<assocMe->pdgId()<<
", energy: "<<assocMe->kineticEnergy()
238 <<
", genParticle: "<<assocMe->genParticleLink().cptr()<<
", beta: "<<simHit->beta()
242 associatedHits.
push_back(std::move(link));
244 int nPrecisionHits = nMdt + nMm + nStgcEta;
245 int nPhiLayers = nTgcPhi + nRpcPhi + nStgcPhi;
247 if (nPrecisionHits < 3) {
251 xAOD::MuonSegment* truthSegment =
out.segments.push_back(std::make_unique<xAOD::MuonSegment>());
252 out.ptDecor(*truthSegment) = muonPt(*simHit, globDir);
257 constexpr
float betaLowLimit = 1.e-6;
262 truthSegment->
setPosition(globPos.x(), globPos.y(), globPos.z());
263 truthSegment->
setDirection(globDir.x(), globDir.y(), globDir.z());
266 truthSegment->
setNHits(nPrecisionHits, nPhiLayers, nTgcEta + nRpcEta);
268 m_idHelperSvc->chamberIndex(segId),
269 m_idHelperSvc->stationEta(segId),
270 m_idHelperSvc->technologyIndex(segId));
272 if (nPhiLayers == 0){
275 truthSegment->
setFitQuality(
chi2, (nPrecisionHits + nPhiLayers + nTgcEta + nRpcEta - 5));
277 out.hitLinkDecor(*truthSegment) = std::move(associatedHits);
285 using HitsPerParticle = std::unordered_map<HepMC::ConstGenParticlePtr, SimHitVec_t>;
286 using HitCollector = std::unordered_map<const MuonGMR4::SpectrometerSector*, HitsPerParticle>;
287 HitCollector hitCollector{};
295 auto genLink = simHit->genParticleLink();
297 if (genLink.isValid()){
298 genParticle = genLink.cptr();
300 if ( (!genParticle && (!m_includePileUpHits || simHit->kineticEnergy() < m_pileUpHitMinE
302 ATH_MSG_VERBOSE(
"Skip hit "<<m_idHelperSvc->toString(simHit->identify())<<
303 " pdgId: "<<simHit->pdgId()<<
", energy: "<<simHit->kineticEnergy()
304 <<
", genParticle: "<<genParticle);
308 hitCollector[
id][genParticle].emplace_back(simHit,
309 toChTrf *xAOD::toEigen(simHit->localPosition()),
310 toChTrf.linear()* xAOD::toEigen(simHit->localDirection()));
315 ATH_CHECK(writeHandle.record(std::make_unique<xAOD::MuonSegmentContainer>(),
316 std::make_unique<xAOD::MuonSegmentAuxContainer>()));
320 for (
auto& [
chamber, collectedParts] : hitCollector) {
322 for (
auto& [
particle, simHits]: collectedParts) {
325 return std::get<1>(
a).z() < std::get<1>(
b).z();
328 buildSegmentsFromBkg(ctx, locToGlob, simHits, writerHolder);
331 constructSegmentFromHits(ctx, locToGlob, simHits, writerHolder);
334 ATH_MSG_DEBUG(
"Constructed "<<writeHandle->size()<<
" truth segments in total ");
335 return StatusCode::SUCCESS;