37 #include "GaudiKernel/IPartPropSvc.h"
45 const float FPGATrackSim_PT_TRUTHMIN = 400.;
46 const float FPGATrackSim_Z_TRUTHMIN = 2300.;
61 SmartIF<IPartPropSvc> partPropSvc{service(
"PartPropSvc")};
83 return StatusCode::SUCCESS;
88 return StatusCode::SUCCESS;
102 event_info.
setLB(eventInfo->lumiBlock());
103 event_info.
setBCID(eventInfo->bcid());
118 std::vector <FPGATrackSimCluster>
clusters;
125 std::vector <FPGATrackSimTruthTrack> truth;
130 std::vector <FPGATrackSimOfflineTrack>
offline;
139 return StatusCode::SUCCESS;
146 ATH_MSG_DEBUG(
"read Offline tracks, size= " << offlineTracksHandle->size());
152 tmpOfflineTrack.
setQOverPt(trackParticle->pt() > 0 ? trackParticle->charge() / trackParticle->pt() : 0);
153 tmpOfflineTrack.
setEta(trackParticle->eta());
154 tmpOfflineTrack.
setPhi(trackParticle->phi());
155 tmpOfflineTrack.
setD0(trackParticle->d0());
156 tmpOfflineTrack.
setZ0(trackParticle->z0());
158 const Trk::TrackStates* trackStates = trackParticle->track()->trackStateOnSurfaces();
159 if (trackStates ==
nullptr) {
161 return StatusCode::FAILURE;
164 if (tsos ==
nullptr)
continue;
167 if (tsos->trackParameters() !=
nullptr &&
168 tsos->trackParameters()->associatedSurface().associatedDetectorElement() !=
nullptr &&
169 tsos->trackParameters()->associatedSurface().associatedDetectorElement()->identify() != 0
193 tmpOfflineHit.
setLocY(-99999.9);
195 tmpOfflineTrack.
addHit(tmpOfflineHit);
199 offline.push_back(tmpOfflineTrack);
203 return StatusCode::SUCCESS;
213 unsigned int hitIndex = 0
u;
218 return StatusCode::SUCCESS;
231 if (pixel_rdoCollection ==
nullptr) {
continue; }
244 hitIndexMap[rdoId] = hitIndex;
249 hitIndexMap[tmpId] = hitIndex;
255 InDetSimDataCollection::const_iterator iter(pixelSDOHandle->find(rdoId));
256 if (
nCells > 1 && iter == pixelSDOHandle->end()) {
258 for (
int ii = 0; ii <
nCells && iter == pixelSDOHandle->end(); ++ii) {
263 if (iter != pixelSDOHandle->end()) { bestTruthLink =
getTruthInformation(iter, parentMask); }
278 else if (barrel_ec == 2)
280 else if (barrel_ec == -2)
293 tmpSGhit.
setToT(pixelRawData->getToT());
297 tmpSGhit.
setBarcode(bestTruthLink->barcode());
306 tmpSGhit.
setBarcodePt(
static_cast<unsigned long>(std::ceil(bestParent ? bestParent->momentum().perp() : 0.)));
319 return StatusCode::SUCCESS;
325 constexpr
int MaxChannelinStripRow = 128;
331 if (SCT_Collection ==
nullptr) {
continue; }
333 std::map<int, bool> firedStrips;
339 const Identifier rdoId = sctRawData->identify();
341 for(
int i = 0;
i < sctRawData->getGroupSize();
i++) {
342 firedStrips[baseLineStrip+
i] =
true;
348 std::map<int, int> stripEncodingForITK;
349 for(
auto& [stripID, fired]: firedStrips)
357 std::bitset<3> hitMap;
361 int currChipID = stripID / MaxChannelinStripRow;
363 int maxStripIDForCurrChip = (currChipID + 1) * MaxChannelinStripRow;
365 for(
int i = 0;
i < 3;
i++)
368 if((stripID + 1 +
i) >= maxStripIDForCurrChip)
continue;
370 if(firedStrips.find(stripID + 1 +
i) != firedStrips.end())
372 if(firedStrips.at(stripID + 1 +
i))
375 firedStrips[stripID + 1 +
i] =
false;
385 stripEncodingForITK[stripID] = (
int)(hitMap.to_ulong());
392 const Identifier rdoId = sctRawData->identify();
396 std::pair<Amg::Vector3D, Amg::Vector3D> endsOfStrip = sielement->
endsOfStrip(LocalPos);
398 hitIndexMap[rdoId] = hitIndex;
405 InDetSimDataCollection::const_iterator iter(stripSDOHandle->find(rdoId));
407 if (iter != stripSDOHandle->end()) { bestTruthLink =
getTruthInformation(iter, parentMask); }
421 else if (barrel_ec == 2)
423 else if (barrel_ec == -2)
436 tmpSGhit.
setBarcode(bestTruthLink->barcode());
447 if(stripEncodingForITK.find(stripID) != stripEncodingForITK.end())
451 int chipID = stripID / MaxChannelinStripRow;
452 int ITkStripID = stripID % MaxChannelinStripRow;
461 ITkStripID +=
offset * MaxChannelinStripRow;
469 tmpSGhit.
setBarcodePt(
static_cast<unsigned long>(std::ceil(bestParent ? bestParent->momentum().perp() : 0.)));
471 tmpSGhit.
setX(0.5 * (endsOfStrip.first.x() + endsOfStrip.second.x()));
472 tmpSGhit.
setY(0.5 * (endsOfStrip.first.y() + endsOfStrip.second.y()));
473 tmpSGhit.
setZ(0.5 * (endsOfStrip.first.z() + endsOfStrip.second.z()));
486 return StatusCode::SUCCESS;
492 unsigned int pixelClusterIndex = 0;
497 if (pixelClusterCollection ==
nullptr) {
507 for (
const Identifier& rdoId : cluster->rdoList()) {
513 InDetSimDataCollection::const_iterator iter(pixelSDOHandle->find(rdoId));
515 if (
nCells > 1 && iter == pixelSDOHandle->end()) {
517 for (
int ii = 0; ii <
nCells && iter == pixelSDOHandle->end(); ++ii) {
525 pixelClusterIndexMap[theId] = pixelClusterIndex;
530 return StatusCode::SUCCESS;
543 if (pixelClusterCollection ==
nullptr) {
547 const int size = pixelClusterCollection->size();
555 for (
const Identifier& rdoId : cluster->rdoList()) {
560 InDetSimDataCollection::const_iterator iter(pixelSDOHandle->find(rdoId));
562 if (
nCells > 1 && iter == pixelSDOHandle->end()) {
564 for (
int ii = 0; ii <
nCells && iter == pixelSDOHandle->end(); ++ii) {
569 if (iter != pixelSDOHandle->end()) { bestTruthLink =
getTruthInformation(iter, parentMask); }
585 clusterEquiv.
setX(globalPos.x());
586 clusterEquiv.
setY(globalPos.y());
587 clusterEquiv.
setZ(globalPos.z());
595 else if (barrel_ec == 2)
597 else if (barrel_ec == -2)
606 clusterEquiv.
setPhiWidth(cluster->width().colRow()[1]);
607 clusterEquiv.
setEtaWidth(cluster->width().colRow()[0]);
611 clusterEquiv.
setBarcode(bestTruthLink->barcode());
620 clusterEquiv.
setBarcodePt(
static_cast<unsigned long>(std::ceil(bestParent ? bestParent->momentum().perp() : 0.)));
635 if (SCT_Collection ==
nullptr) {
continue; }
637 const Identifier rdoId = sctRawData->identify();
648 InDetSimDataCollection::const_iterator iter(stripSDOHandle->find(rdoId));
650 if (iter != stripSDOHandle->end()) { bestTruthLink =
getTruthInformation(iter, parentMask); }
668 else if (barrel_ec == 2)
670 else if (barrel_ec == -2)
680 clusterEquiv.
setPhiWidth(sctRawData->getGroupSize());
684 clusterEquiv.
setBarcode(bestTruthLink->barcode());
692 clusterEquiv.
setBarcodePt(
static_cast<unsigned long>(std::ceil(bestParent ? bestParent->momentum().perp() : 0.)));
700 return StatusCode::SUCCESS;
707 ATH_MSG_DEBUG(
"Dump truth tracks, size " << simTracksHandle->size());
710 for (
unsigned int ievt = 0; ievt < simTracksHandle->size(); ++ievt) {
711 const HepMC::GenEvent* genEvent = simTracksHandle->at(ievt);
713 HepGeom::Point3D<double> primaryVtx(0., 0., 0.);
718 primaryVtx.set(spv->position().x(),
720 spv->position().z());
721 ATH_MSG_DEBUG(
"using signal process vertex for eventIndex " << ievt <<
":"
722 << primaryVtx.x() <<
"\t" << primaryVtx.y() <<
"\t" << primaryVtx.z());
724 for (
const auto&
particle: *genEvent) {
725 const int pdgcode =
particle->pdg_id();
727 if (
particle->production_vertex() ==
nullptr) {
736 if (pdgcode < 0)
charge *= -1.;
737 if (std::abs(
charge) < 0.5) {
757 const double track_truth_d0 = tP ? tP->parameters()[
Trk::d0] : 999.;
758 const double track_truth_phi = tP ? tP->parameters()[
Trk::phi] : 999.;
759 const double track_truth_p = (tP && fabs(tP->parameters()[
Trk::qOverP]) > 1.e-8) ?
761 const double track_truth_x0 = tP ? tP->
position().x() : 999.;
762 const double track_truth_y0 = tP ? tP->
position().y() : 999.;
763 const double track_truth_z0 = tP ? tP->parameters()[
Trk::z0] : 999.;
764 const double track_truth_q = tP ? tP->
charge() : 0.;
765 const double track_truth_sinphi = tP ?
std::sin(tP->parameters()[
Trk::phi]) : -1.;
766 const double track_truth_cosphi = tP ?
std::cos(tP->parameters()[
Trk::phi]) : -1.;
767 const double track_truth_sintheta = tP ?
std::sin(tP->parameters()[
Trk::theta]) : -1.;
768 const double track_truth_costheta = tP ?
std::cos(tP->parameters()[
Trk::theta]) : -1.;
769 double truth_d0corr = track_truth_d0 - (primaryVtx.y() *
cos(track_truth_phi) - primaryVtx.x() *
sin(track_truth_phi));
770 double truth_zvertex = 0.;
771 const HepGeom::Point3D<double> startVertex(
particle->production_vertex()->position().x(),
particle->production_vertex()->position().y(),
particle->production_vertex()->position().z());
774 if (std::abs(truth_d0corr) > 2.) {
isPrimary =
false; }
779 const HepGeom::Point3D<double> startVertex(
particle->production_vertex()->position().x(),
particle->production_vertex()->position().y(),
particle->production_vertex()->position().z());
780 if (std::abs(startVertex.z() - truth_zvertex) > 100.) {
isPrimary =
false; }
782 HepGeom::Point3D<double> endVertex(
particle->end_vertex()->position().x(),
particle->end_vertex()->position().y(),
particle->end_vertex()->position().z());
783 if (endVertex.perp() < FPGATrackSim_PT_TRUTHMIN && std::abs(endVertex.z()) < FPGATrackSim_Z_TRUTHMIN) {
isPrimary =
false; }
793 tmpSGTrack.
setVtxX(track_truth_x0);
794 tmpSGTrack.
setVtxY(track_truth_y0);
795 tmpSGTrack.
setVtxZ(track_truth_z0);
796 tmpSGTrack.
setD0(track_truth_d0);
797 tmpSGTrack.
setZ0(track_truth_z0);
798 tmpSGTrack.
setVtxZ(primaryVtx.z());
799 tmpSGTrack.
setQ(track_truth_q);
800 tmpSGTrack.
setPX(track_truth_p * (track_truth_cosphi * track_truth_sintheta));
801 tmpSGTrack.
setPY(track_truth_p * (track_truth_sinphi * track_truth_sintheta));
802 tmpSGTrack.
setPZ(track_truth_p * track_truth_costheta);
810 truth.push_back(tmpSGTrack);
815 return StatusCode::SUCCESS;
823 const std::vector<InDetSimData::Deposit>& deposits(sdo.
getdeposits());
824 float bestPt{-999.f};
829 if (!particleLink.
isValid()) {
continue; }
830 const float genEta = particleLink->momentum().pseudoRapidity();
831 const float genPt = particleLink->momentum().perp();
837 if (std::fabs(genEta) >
m_maxEta) {
continue; }
839 if (bestPt < genPt) {
841 bestTruthLink = &particleLink;
850 return bestTruthLink;