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)
295 tmpSGhit.
setToT(pixelRawData->getToT());
299 tmpSGhit.
setBarcode(bestTruthLink->barcode());
308 tmpSGhit.
setBarcodePt(
static_cast<unsigned long>(std::ceil(bestParent ? bestParent->momentum().perp() : 0.)));
321 return StatusCode::SUCCESS;
327 constexpr
int MaxChannelinStripRow = 128;
333 if (SCT_Collection ==
nullptr) {
continue; }
335 std::map<int, bool> firedStrips;
341 const Identifier rdoId = sctRawData->identify();
343 for(
int i = 0;
i < sctRawData->getGroupSize();
i++) {
344 firedStrips[baseLineStrip+
i] =
true;
350 std::map<int, int> stripEncodingForITK;
351 for(
auto& [stripID, fired]: firedStrips)
359 std::bitset<3> hitMap;
363 int currChipID = stripID / MaxChannelinStripRow;
365 int maxStripIDForCurrChip = (currChipID + 1) * MaxChannelinStripRow;
367 for(
int i = 0;
i < 3;
i++)
370 if((stripID + 1 +
i) >= maxStripIDForCurrChip)
continue;
372 if(firedStrips.find(stripID + 1 +
i) != firedStrips.end())
374 if(firedStrips.at(stripID + 1 +
i))
377 firedStrips[stripID + 1 +
i] =
false;
387 stripEncodingForITK[stripID] = (
int)(hitMap.to_ulong());
394 const Identifier rdoId = sctRawData->identify();
398 std::pair<Amg::Vector3D, Amg::Vector3D> endsOfStrip = sielement->
endsOfStrip(localPos);
400 hitIndexMap[rdoId] = hitIndex;
407 InDetSimDataCollection::const_iterator iter(stripSDOHandle->find(rdoId));
409 if (iter != stripSDOHandle->end()) { bestTruthLink =
getTruthInformation(iter, parentMask); }
423 else if (barrel_ec == 2)
425 else if (barrel_ec == -2)
440 tmpSGhit.
setBarcode(bestTruthLink->barcode());
451 if(stripEncodingForITK.find(stripID) != stripEncodingForITK.end())
455 int chipID = stripID / MaxChannelinStripRow;
456 int ITkStripID = stripID % MaxChannelinStripRow;
465 ITkStripID +=
offset * MaxChannelinStripRow;
473 tmpSGhit.
setBarcodePt(
static_cast<unsigned long>(std::ceil(bestParent ? bestParent->momentum().perp() : 0.)));
475 tmpSGhit.
setX(0.5 * (endsOfStrip.first.x() + endsOfStrip.second.x()));
476 tmpSGhit.
setY(0.5 * (endsOfStrip.first.y() + endsOfStrip.second.y()));
477 tmpSGhit.
setZ(0.5 * (endsOfStrip.first.z() + endsOfStrip.second.z()));
490 return StatusCode::SUCCESS;
496 unsigned int pixelClusterIndex = 0;
501 if (pixelClusterCollection ==
nullptr) {
511 for (
const Identifier& rdoId : cluster->rdoList()) {
517 InDetSimDataCollection::const_iterator iter(pixelSDOHandle->find(rdoId));
519 if (
nCells > 1 && iter == pixelSDOHandle->end()) {
521 for (
int ii = 0; ii <
nCells && iter == pixelSDOHandle->end(); ++ii) {
529 pixelClusterIndexMap[theId] = pixelClusterIndex;
534 return StatusCode::SUCCESS;
547 if (pixelClusterCollection ==
nullptr) {
551 const int size = pixelClusterCollection->size();
559 for (
const Identifier& rdoId : cluster->rdoList()) {
564 InDetSimDataCollection::const_iterator iter(pixelSDOHandle->find(rdoId));
566 if (
nCells > 1 && iter == pixelSDOHandle->end()) {
568 for (
int ii = 0; ii <
nCells && iter == pixelSDOHandle->end(); ++ii) {
573 if (iter != pixelSDOHandle->end()) { bestTruthLink =
getTruthInformation(iter, parentMask); }
589 clusterEquiv.
setX(globalPos.x());
590 clusterEquiv.
setY(globalPos.y());
591 clusterEquiv.
setZ(globalPos.z());
599 else if (barrel_ec == 2)
601 else if (barrel_ec == -2)
612 clusterEquiv.
setPhiWidth(cluster->width().colRow()[1]);
613 clusterEquiv.
setEtaWidth(cluster->width().colRow()[0]);
617 clusterEquiv.
setBarcode(bestTruthLink->barcode());
626 clusterEquiv.
setBarcodePt(
static_cast<unsigned long>(std::ceil(bestParent ? bestParent->momentum().perp() : 0.)));
641 if (SCT_Collection ==
nullptr) {
continue; }
643 const Identifier rdoId = sctRawData->identify();
654 InDetSimDataCollection::const_iterator iter(stripSDOHandle->find(rdoId));
656 if (iter != stripSDOHandle->end()) { bestTruthLink =
getTruthInformation(iter, parentMask); }
674 else if (barrel_ec == 2)
676 else if (barrel_ec == -2)
688 clusterEquiv.
setPhiWidth(sctRawData->getGroupSize());
692 clusterEquiv.
setBarcode(bestTruthLink->barcode());
701 clusterEquiv.
setBarcodePt(
static_cast<unsigned long>(std::ceil(bestParent ? bestParent->momentum().perp() : 0.)));
709 return StatusCode::SUCCESS;
716 ATH_MSG_DEBUG(
"Dump truth tracks, size " << simTracksHandle->size());
719 for (
unsigned int ievt = 0; ievt < simTracksHandle->size(); ++ievt) {
720 const HepMC::GenEvent* genEvent = simTracksHandle->at(ievt);
722 HepGeom::Point3D<double> primaryVtx(0., 0., 0.);
727 primaryVtx.set(spv->position().x(),
729 spv->position().z());
730 ATH_MSG_DEBUG(
"using signal process vertex for eventIndex " << ievt <<
":"
731 << primaryVtx.x() <<
"\t" << primaryVtx.y() <<
"\t" << primaryVtx.z());
733 for (
const auto&
particle: *genEvent) {
734 const int pdgcode =
particle->pdg_id();
736 if (
particle->production_vertex() ==
nullptr) {
745 if (pdgcode < 0)
charge *= -1.;
746 if (std::abs(
charge) < 0.5) {
766 const double track_truth_d0 = tP ? tP->parameters()[
Trk::d0] : 999.;
767 const double track_truth_phi = tP ? tP->parameters()[
Trk::phi] : 999.;
768 const double track_truth_p = (tP && fabs(tP->parameters()[
Trk::qOverP]) > 1.e-8) ?
770 const double track_truth_x0 = tP ? tP->
position().x() : 999.;
771 const double track_truth_y0 = tP ? tP->
position().y() : 999.;
772 const double track_truth_z0 = tP ? tP->parameters()[
Trk::z0] : 999.;
773 const double track_truth_q = tP ? tP->
charge() : 0.;
774 const double track_truth_sinphi = tP ?
std::sin(tP->parameters()[
Trk::phi]) : -1.;
775 const double track_truth_cosphi = tP ?
std::cos(tP->parameters()[
Trk::phi]) : -1.;
776 const double track_truth_sintheta = tP ?
std::sin(tP->parameters()[
Trk::theta]) : -1.;
777 const double track_truth_costheta = tP ?
std::cos(tP->parameters()[
Trk::theta]) : -1.;
778 double truth_d0corr = track_truth_d0 - (primaryVtx.y() *
cos(track_truth_phi) - primaryVtx.x() *
sin(track_truth_phi));
779 double truth_zvertex = 0.;
780 const HepGeom::Point3D<double> startVertex(
particle->production_vertex()->position().x(),
particle->production_vertex()->position().y(),
particle->production_vertex()->position().z());
783 if (std::abs(truth_d0corr) > 2.) {
isPrimary =
false; }
788 const HepGeom::Point3D<double> startVertex(
particle->production_vertex()->position().x(),
particle->production_vertex()->position().y(),
particle->production_vertex()->position().z());
789 if (std::abs(startVertex.z() - truth_zvertex) > 100.) {
isPrimary =
false; }
791 HepGeom::Point3D<double> endVertex(
particle->end_vertex()->position().x(),
particle->end_vertex()->position().y(),
particle->end_vertex()->position().z());
792 if (endVertex.perp() < FPGATrackSim_PT_TRUTHMIN && std::abs(endVertex.z()) < FPGATrackSim_Z_TRUTHMIN) {
isPrimary =
false; }
802 tmpSGTrack.
setVtxX(track_truth_x0);
803 tmpSGTrack.
setVtxY(track_truth_y0);
804 tmpSGTrack.
setVtxZ(track_truth_z0);
805 tmpSGTrack.
setD0(track_truth_d0);
806 tmpSGTrack.
setZ0(track_truth_z0);
807 tmpSGTrack.
setVtxZ(primaryVtx.z());
808 tmpSGTrack.
setQ(track_truth_q);
809 tmpSGTrack.
setPX(track_truth_p * (track_truth_cosphi * track_truth_sintheta));
810 tmpSGTrack.
setPY(track_truth_p * (track_truth_sinphi * track_truth_sintheta));
811 tmpSGTrack.
setPZ(track_truth_p * track_truth_costheta);
819 truth.push_back(tmpSGTrack);
824 return StatusCode::SUCCESS;
832 const std::vector<InDetSimData::Deposit>& deposits(sdo.
getdeposits());
833 float bestPt{-999.f};
838 if (!particleLink.
isValid()) {
continue; }
839 const float genEta = particleLink->momentum().pseudoRapidity();
840 const float genPt = particleLink->momentum().perp();
846 if (std::fabs(genEta) >
m_maxEta) {
continue; }
848 if (bestPt < genPt) {
850 bestTruthLink = &particleLink;
859 return bestTruthLink;