29 const std::vector<float> emptyVec;
32 const std::set<int> bad_origins{
42 constexpr
float dummy_val = -1.;
45 RecordPars() =
default;
48 mom{std::move(_mom)}{}
49 RecordPars(
const CLHEP::Hep3Vector& _pos,
const CLHEP::Hep3Vector& _mom):
50 pos{_pos.x(), _pos.y(), _pos.z()},
51 mom{_mom.x(), _mom.y(), _mom.z()}{
57 std::string record_name;
81 return StatusCode::SUCCESS;
88 if (!truthContainer.
isPresent())
return StatusCode::SUCCESS;
89 if (!truthContainer.
isValid()) {
91 return StatusCode::FAILURE;
96 ATH_CHECK(muonTruthContainer.
record(std::make_unique<xAOD::TruthParticleContainer>(),
97 std::make_unique<xAOD::TruthParticleAuxContainer>()));
103 segmentContainer.
record(std::make_unique<xAOD::MuonSegmentContainer>(), std::make_unique<xAOD::MuonSegmentAuxContainer>()));
104 ATH_MSG_DEBUG(
"Recorded MuonSegmentContainer with key: " << segmentContainer.
name());
109 if (!
MC::isStable(truth) || !truth->isMuon() || truth->pt() < 1000.)
continue;
111 muonTruthContainer->
push_back(truthParticle);
112 truthParticle->
setPdgId(truth->pdgId());
114 truthParticle->
setStatus(truth->status());
115 truthParticle->
setPx(truth->px());
116 truthParticle->
setPy(truth->py());
117 truthParticle->
setPz(truth->pz());
118 truthParticle->
setE(truth->e());
119 truthParticle->
setM(truth->m());
120 if (truth->hasProdVtx()) truthParticle->
setProdVtxLink(truth->prodVtxLink());
123 ATH_MSG_DEBUG(
"Found stable muon: " << truth->pt() <<
" eta " << truth->eta() <<
" phi " << truth->phi() <<
" mass "
124 << truth->m() <<
" barcode " <<
HepMC::barcode(truth) <<
" truthParticle->barcode "
133 std::pair<MCTruthPartClassifier::ParticleType, MCTruthPartClassifier::ParticleOrigin> truthClass =
135 iType = truthClass.first;
136 iOrigin = truthClass.second;
138 dec_truthOrigin(*truthParticle) = iOrigin;
139 dec_truthType(*truthParticle) = iType;
152 bool goodMuon = bad_origins.find(iOrigin) == bad_origins.end();
154 ATH_MSG_DEBUG(
"good muon with type " << iType <<
" and origin" << iOrigin);
165 return StatusCode::SUCCESS;
174 std::vector<SG::ReadHandle<MuonSimDataCollection> > sdoCollections(6);
177 if (!
col.isPresent()) {
178 ATH_MSG_ERROR(
"MuonSimDataCollection " <<
col.name() <<
" not in StoreGate");
179 return StatusCode::FAILURE;
186 if (
index >= (
int)sdoCollections.size()) {
189 sdoCollections[
index] = std::move(
col);
195 std::map<Muon::MuonStationIndex::ChIndex, int> matchMap;
198 for (
const auto& lay :
ids) {
201 bool firstPosSet{
false}, secondPosSet{
false};
204 uint8_t nprecLayers{0}, nphiLayers{0}, ntrigEtaLayers{0};
205 std::set<int> phiLayers, etaLayers, precLayers;
208 for (
const auto&
id : lay.second) {
224 precLayers.insert(iid);
239 if (
index < (
int)sdoCollections.size() && !sdoCollections[
index]->empty()) {
240 auto pos = sdoCollections[
index]->find(
id);
242 gpos =
pos->second.globalPosition();
243 if (gpos.perp() > 0.1) ok =
true;
252 return std::abs(
p1.z()) < std::abs(
p2.z());
254 return p1.perp() <
p2.perp();
259 }
else if (!secondPosSet) {
262 if (isSmaller(secondPos, firstPos))
std::swap(firstPos, secondPos);
265 if (isSmaller(gpos, firstPos))
267 else if (isSmaller(secondPos, gpos))
272 auto pos = cscCollection->find(
id);
273 if (
pos == cscCollection->end()) {
277 ATH_MSG_DEBUG(
"found csc sdo with " <<
pos->second.getdeposits().size() <<
" deposits");
278 Amg::Vector3D locpos(0,
pos->second.getdeposits()[0].second.ypos(),
pos->second.getdeposits()[0].second.zpos());
284 }
else if (!secondPosSet) {
287 if (secondPos.perp() < firstPos.perp())
std::swap(firstPos, secondPos);
289 if (gpos.perp() < firstPos.perp())
291 else if (secondPos.perp() < gpos.perp())
296 if (precLayers.size() > 2) {
297 matchMap[lay.first] =
index;
298 if (!phiLayers.empty()) nphiLayers = phiLayers.size();
299 ntrigEtaLayers = etaLayers.size();
300 nprecLayers = precLayers.size();
301 ATH_MSG_DEBUG(
" total counts: precision " <<
static_cast<int>(nprecLayers) <<
" phi layers " <<
static_cast<int>(nphiLayers)
302 <<
" eta trig layers " <<
static_cast<int>(ntrigEtaLayers)
304 <<
" truthLink " << truthLink);
307 segment->setNHits(nprecLayers, nphiLayers, ntrigEtaLayers);
309 truthParticleLinkAcc(
"truthParticleLink");
310 truthParticleLinkAcc(*
segment) = truthLink;
316 segment->setIdentifier(sector, chIndex, eta, technology);
318 if (firstPosSet && secondPosSet) {
321 ATH_MSG_DEBUG(
" got position : r " << gpos.perp() <<
" z " << gpos.z() <<
" and direction: theta " << gdir.theta()
322 <<
" phi " << gdir.phi());
323 segment->setPosition(gpos.x(), gpos.y(), gpos.z());
324 segment->setDirection(gdir.x(), gdir.y(), gdir.z());
328 return StatusCode::SUCCESS;
342 if (!trackRecordCollection.isPresent()) {
343 ATH_MSG_FATAL(
"Failed to retrieve "<<trackRecordCollection.key());
344 return StatusCode::FAILURE;
346 const std::string r_name = trackRecordCollection.key();
354 float&
x = xAcc(truthParticle);
355 float&
y = yAcc(truthParticle);
356 float&
z = zAcc(truthParticle);
357 float&
px = pxAcc(truthParticle);
358 float&
py = pyAcc(truthParticle);
359 float&
pz = pzAcc(truthParticle);
362 bool& found_truth = matchedAcc(truthParticle);
374 float& ex = exAcc(truthParticle);
375 float& ey = eyAcc(truthParticle);
376 float& ez = ezAcc(truthParticle);
377 float& epx = epxAcc(truthParticle);
378 float& epy = epyAcc(truthParticle);
379 float& epz = epzAcc(truthParticle);
382 ecovAcc(truthParticle) = emptyVec;
385 eisAcc(truthParticle) =
false;
386 ex = ey = ez = epx = epy = epz = dummy_val;
389 for (
const auto& trackRecord : *trackRecordCollection) {
391 CLHEP::Hep3Vector
pos = trackRecord.GetPosition();
392 CLHEP::Hep3Vector
mom = trackRecord.GetMomentum();
393 ATH_MSG_VERBOSE(
"Found associated " << r_name <<
" pt " <<
mom.perp() <<
" position: r " <<
pos.perp() <<
" z " <<
pos.z());
401 if (!trackingGeometry)
continue;
406 r_pars.record_name = r_name;
408 if (r_name ==
"CaloEntryLayer")
409 volume = trackingGeometry->
trackingVolume(
"InDet::Containers::InnerDetector");
410 else if (r_name ==
"MuonEntryLayer")
412 else if (r_name ==
"MuonExitLayer")
413 volume = trackingGeometry->
trackingVolume(
"Muon::Containers::MuonSystem");
417 r_pars.volume = found_truth ? volume :
nullptr;
422 if (!trackingGeometry)
return StatusCode::SUCCESS;
428 cov(4, 4) = 1
e-3 / truthParticle.
p4().P();
433 if ( (!start_pars.record_name.empty() && !start_pars.volume) || !end_pars.volume)
continue;
437 const std::string& r_name = end_pars.record_name;
444 float& ex = exAcc(truthParticle);
445 float& ey = eyAcc(truthParticle);
446 float& ez = ezAcc(truthParticle);
447 float& epx = epxAcc(truthParticle);
448 float& epy = epyAcc(truthParticle);
449 float& epz = epzAcc(truthParticle);
452 std::vector<float>& covMat = ecovAcc(truthParticle);
454 std::unique_ptr<Trk::TrackParameters> exPars{
461 eisAcc(truthParticle) =
true;
462 ex = exPars->position().x();
463 ey = exPars->position().y();
464 ez = exPars->position().z();
465 epx = exPars->momentum().x();
466 epy = exPars->momentum().y();
467 epz = exPars->momentum().z();
472 <<
" truth: r " << end_pars.pos.perp() <<
" z "
473 << end_pars.pos.z() <<
" p "
474 << end_pars.mom.mag() << std::endl
475 <<
" extrp: r " << exPars->position().perp() <<
" z "
476 << exPars->position().z() <<
" p "
477 << exPars->momentum().mag() <<
" res p "
478 << (end_pars.mom.mag() -
479 exPars->momentum().mag())
480 <<
" error " << errorp <<
" cov "
483 << (end_pars.mom.mag() -
484 exPars->momentum().mag()) /
487 return StatusCode::SUCCESS;
494 std::vector<unsigned int> nprecHitsPerChamberLayer;
496 std::vector<unsigned int> nphiHitsPerChamberLayer;
498 std::vector<unsigned int> ntrigEtaHitsPerChamberLayer;
504 if (!
col.isPresent()) {
505 ATH_MSG_FATAL(
"PRD_MultiTruthCollection " <<
col.name() <<
" not in StoreGate");
506 return StatusCode::FAILURE;
509 for (
const auto& trajectory : *
col) {
511 if (
std::find(truthParticleHistory.begin(),truthParticleHistory.end(),
HepMC::barcode(trajectory.second)) == truthParticleHistory.end())
continue;
528 ids[chIndex].push_back(
id);
534 ++nphiHitsPerChamberLayer.at(
index);
536 ++nprecHitsPerChamberLayer.at(chIndex);
539 ++nprecHitsPerChamberLayer.at(chIndex);
544 ++nphiHitsPerChamberLayer.at(
index);
546 ++ntrigEtaHitsPerChamberLayer.at(
index);
551 ++nphiHitsPerChamberLayer.at(
index);
554 ++nprecHitsPerChamberLayer.at(chIndex);
672 nprecLayersAcc(truthParticle) = nprecLayers;
673 nphiLayersAcc(truthParticle) = nphiLayers;
674 ntrigEtaLayersAcc(truthParticle) = ntrigEtaLayers;
703 ATH_MSG_DEBUG(
"Precision layers " <<
static_cast<int>(nprecLayers) <<
" phi layers " <<
static_cast<int>(nphiLayers)
704 <<
" triggerEta layers " <<
static_cast<int>(ntrigEtaLayers));
706 if (nprecLayers > 0) {
709 for (
int index = 0; index < static_cast<int>(nprecHitsPerChamberLayer.size()); ++
index) {
710 if (nprecHitsPerChamberLayer[
index] > 0)
712 <<
" hits " << nprecHitsPerChamberLayer[
index];
715 if (nphiLayers > 0) {
717 for (
int index = 0; index < static_cast<int>(nphiHitsPerChamberLayer.size()); ++
index) {
718 if (nphiHitsPerChamberLayer[
index] > 0)
720 <<
" hits " << nphiHitsPerChamberLayer[
index];
724 if (ntrigEtaLayers > 0) {
726 for (
int index = 0; index < static_cast<int>(ntrigEtaHitsPerChamberLayer.size()); ++
index) {
727 if (ntrigEtaHitsPerChamberLayer[
index] > 0)
729 <<
" hits " << ntrigEtaHitsPerChamberLayer[
index];
734 return StatusCode::SUCCESS;
741 truthMdtHitsAcc(
"truthMdtHits");
743 truthTgcHitsAcc(
"truthTgcHits");
745 truthRpcHitsAcc(
"truthRpcHits");
746 std::vector<unsigned long long>& mdtTruthHits = truthMdtHitsAcc(truthParticle);
747 std::vector<unsigned long long>& tgcTruthHits = truthTgcHitsAcc(truthParticle);
748 std::vector<unsigned long long>& rpcTruthHits = truthRpcHitsAcc(truthParticle);
750 std::vector<unsigned long long> stgcTruthHits;
751 std::vector<unsigned long long> cscTruthHits;
752 std::vector<unsigned long long> mmTruthHits;
755 int nEI = 0, nEM = 0;
756 for (
const auto& lay :
ids) {
760 for (
const auto&
id : lay.second) {
762 mdtTruthHits.push_back(
id.get_compact());
764 cscTruthHits.push_back(
id.get_compact());
770 tgcTruthHits.push_back(
id.get_compact());
772 stgcTruthHits.push_back(
id.get_compact());
774 rpcTruthHits.push_back(
id.get_compact());
776 mmTruthHits.push_back(
id.get_compact());
781 truthCscHitsAcc(
"truthCscHits");
782 truthCscHitsAcc(truthParticle) = std::move(cscTruthHits);
786 truthStgcHitsAcc(
"truthStgcHits");
787 truthStgcHitsAcc(truthParticle) = std::move(stgcTruthHits);
791 truthMMHitsAcc(
"truthMMHits");
792 truthMMHitsAcc(truthParticle) = std::move(mmTruthHits);
794 ATH_MSG_VERBOSE(
"Added " << mdtTruthHits.size() <<
" mdt truth hits, " << cscTruthHits.size() <<
" csc truth hits, "
795 << rpcTruthHits.size() <<
" rpc truth hits, and " << tgcTruthHits.size() <<
" tgc truth hits");