372 const EventContext &ctx = Gaudi::Hive::currentContext();
378 std::map<Identifier, long int> clusterIDMapIdx;
381 std::map<Identifier, long int> clusterIDMapSpacePointIdx;
385 std::map<std::pair<int, int>, std::pair<bool, int>> allTruthParticles;
389 if (not mcEventCollectionHandle.isValid()) {
391 return StatusCode::FAILURE;
393 mcCollptr = mcEventCollectionHandle.cptr();
398 if (not eventInfoHandle.isValid()) {
400 return StatusCode::FAILURE;
402 eventInfo = eventInfoHandle.cptr();
407 std::map<int, int> allSubEvents;
411 bool duplicateSubeventID =
false;
412 for (
unsigned int cntr = 0; cntr < mcCollptr->
size(); ++cntr) {
413 int ID = mcCollptr->
at(cntr)->event_number();
422 if (
it == allSubEvents.end())
423 allSubEvents.insert(std::make_pair(
ID, 1));
426 duplicateSubeventID =
true;
430 if (duplicateSubeventID) {
437 (*m_Part_vParentID).clear();
438 (*m_Part_vParentBarcode).clear();
441 for (
unsigned int cntr = 0; cntr < mcCollptr->
size(); ++cntr) {
442 const HepMC::GenEvent *genEvt = (mcCollptr->
at(cntr));
451 for (
auto p : *genEvt) {
453 float px,
py,
pz,
pt,
eta, vx, vy, vz,
radius,
status,
charge = 0.;
454 std::vector<int> vParentID;
455 std::vector<int> vParentBarcode;
457 int vProdNin, vProdNout, vProdStatus, vProdBarcode;
458 bool passed =
isPassed(
p,
px,
py,
pz,
pt, eta, vx, vy, vz,
radius,
status,
charge, vParentID, vParentBarcode,
459 vProdNin, vProdNout, vProdStatus, vProdBarcode);
460 allTruthParticles.insert(std::make_pair(std::make_pair(genEvt->event_number(),
HepMC::barcode(
p)),
461 std::make_pair(
passed, 0)));
483 (*m_Part_vParentID).push_back(vParentID);
484 (*m_Part_vParentBarcode).push_back(vParentBarcode);
497 if (not pixelClusterContainerHandle.isValid()) {
499 return StatusCode::FAILURE;
505 if (not stripClusterContainerHandle.isValid()) {
507 return StatusCode::FAILURE;
511 auto cartesion_to_spherical = [](
const Amg::Vector3D &xyzVec,
float &eta_,
float &phi_) {
514 r3 += xyzVec[
idx] * xyzVec[
idx];
517 phi_ = atan2(xyzVec[1], xyzVec[0]);
518 float theta_ = acos(xyzVec[2] / r3);
519 eta_ =
log(
tan(0.5 * theta_));
528 (*m_CLhardware).clear();
529 (*m_CLparticleLink_eventIndex).clear();
530 (*m_CLparticleLink_barcode).clear();
531 (*m_CLbarcodesLinked).clear();
532 (*m_CLparticle_charge).clear();
536 (*m_CLlocal_cov).clear();
543 if (not sdoCollectionHandle.isValid()) {
545 return StatusCode::FAILURE;
547 sdoCollection = sdoCollectionHandle.cptr();
551 if (clusterCollection->empty())
562 float norm_x = fabs(my_normal.x()) > 1
e-5 ? my_normal.x() : 0.;
563 float norm_y = fabs(my_normal.y()) > 1
e-5 ? my_normal.y() : 0.;
564 float norm_z = fabs(my_normal.z()) > 1
e-5 ? my_normal.z() : 0.;
569 ATH_MSG_ERROR(
"Dynamic cast failed at " << __LINE__ <<
" of MergedPixelsTool.cxx.");
570 return StatusCode::FAILURE;
574 for (
const auto cluster : *clusterCollection) {
580 const Amg::MatrixX &local_cov = cluster->localCovariance();
582 std::vector<std::pair<int, int>>
barcodes = {};
583 std::vector<int> particleLink_eventIndex = {};
584 std::vector<int> particleLink_barcode = {};
585 std::vector<bool> barcodesLinked = {};
586 std::vector<float>
charge = {};
587 std::vector<int> phis = {};
588 std::vector<int>
etas = {};
589 std::vector<int> tots = {};
595 float charge_count = 0;
598 for (
unsigned int rdo = 0; rdo < cluster->rdoList().
size(); rdo++) {
599 const auto &rdoID = cluster->rdoList().at(rdo);
612 charge_count += cluster->totList().at(rdo);
616 tots.push_back(cluster->totList().at(rdo));
618 auto pos = sdoCollection->find(rdoID);
619 if (
pos != sdoCollection->end()) {
620 for (
auto deposit :
pos->second.getdeposits()) {
626 particleLink_eventIndex.push_back(particleLink.
eventIndex());
627 particleLink_barcode.push_back(particleLink.
barcode());
628 charge.push_back(deposit.second);
629 barcodesLinked.push_back(particleLink.
isValid());
646 Amg::Vector3D localDirection = localEndPosition - localStartPosition;
648 float loc_eta = 0, loc_phi = 0;
649 cartesion_to_spherical(localDirection, loc_eta, loc_phi);
654 Amg::Vector3D direction = globalEndPosition - globalStartPosition;
655 float glob_eta = 0, glob_phi = 0;
656 cartesion_to_spherical(direction, glob_eta, glob_phi);
661 float trkphicomp = direction.dot(my_phiax);
662 float trketacomp = direction.dot(my_etaax);
663 float trknormcomp = direction.dot(my_normal);
664 double phi_angle = atan2(trknormcomp, trkphicomp);
665 double eta_angle = atan2(trknormcomp, trketacomp);
667 clusterIDMapIdx[cluster->identify()] =
m_selected;
668 std::vector<double> v_local_cov;
669 if (local_cov.size() > 0) {
670 for (
size_t i = 0,
nRows = local_cov.rows(),
nCols = local_cov.cols();
i <
nRows;
i++) {
671 for (
size_t j = 0; j <
nCols; ++j) {
672 v_local_cov.push_back(local_cov(
i, j));
679 (*m_CLhardware).push_back(
"PIXEL");
689 (*m_CLparticleLink_eventIndex).push_back(particleLink_eventIndex);
690 (*m_CLparticleLink_barcode).push_back(particleLink_barcode);
691 (*m_CLbarcodesLinked).push_back(barcodesLinked);
692 (*m_CLparticle_charge).push_back(
charge);
693 (*m_CLetas).push_back(
etas);
694 (*m_CLphis).push_back(phis);
695 (*m_CLtots).push_back(tots);
713 (*m_CLlocal_cov).push_back(v_local_cov);
732 if (not sdoCollectionHandle.isValid()) {
734 return StatusCode::FAILURE;
736 sdoCollection = sdoCollectionHandle.cptr();
740 if (clusterCollection->empty())
752 float norm_x = fabs(my_normal.x()) > 1
e-5 ? my_normal.x() : 0.;
753 float norm_y = fabs(my_normal.y()) > 1
e-5 ? my_normal.y() : 0.;
754 float norm_z = fabs(my_normal.z()) > 1
e-5 ? my_normal.z() : 0.;
757 for (
const auto cluster : *clusterCollection) {
763 const Amg::MatrixX &local_cov = cluster->localCovariance();
765 std::vector<std::pair<int, int>>
barcodes = {};
766 std::vector<int> particleLink_eventIndex = {};
767 std::vector<int> particleLink_barcode = {};
768 std::vector<bool> barcodesLinked = {};
769 std::vector<float>
charge = {};
771 std::vector<int> tots = {};
772 std::vector<int> strip_ids = {};
774 int max_strip = -999;
776 float charge_count = 0;
779 for (
unsigned int rdo = 0; rdo < cluster->rdoList().
size(); rdo++) {
780 const auto &rdoID = cluster->rdoList().at(rdo);
784 if (min_strip > strip)
786 if (max_strip < strip)
788 strip_ids.push_back(strip);
793 auto pos = sdoCollection->find(rdoID);
794 if (
pos != sdoCollection->end()) {
795 for (
auto deposit :
pos->second.getdeposits()) {
802 particleLink_eventIndex.push_back(particleLink.
eventIndex());
803 particleLink_barcode.push_back(particleLink.
barcode());
804 charge.push_back(deposit.second);
805 barcodesLinked.push_back(particleLink.
isValid());
815 ATH_MSG_ERROR(
"Failed at " << __LINE__ <<
" of accessing SCT ModuleSide Design");
816 return StatusCode::FAILURE;
820 std::pair<Amg::Vector3D, Amg::Vector3D> ends(
834 Amg::Vector3D localDirection = localEndPosition - localStartPosition;
835 float loc_eta = 0, loc_phi = 0;
836 cartesion_to_spherical(localDirection, loc_eta, loc_phi);
841 Amg::Vector3D direction = globalEndPosition - globalStartPosition;
842 float glob_eta = 0, glob_phi = 0;
843 cartesion_to_spherical(direction, glob_eta, glob_phi);
848 float trkphicomp = direction.dot(my_phiax);
849 float trketacomp = direction.dot(my_etaax);
850 float trknormcomp = direction.dot(my_normal);
851 double phi_angle = atan2(trknormcomp, trkphicomp);
852 double eta_angle = atan2(trknormcomp, trketacomp);
855 clusterIDMapIdx[cluster->identify()] =
m_selected;
857 std::vector<int> cst;
858 for (
unsigned strip = 0;
strip < strip_ids.size();
strip++) {
861 std::vector<double> v_local_cov;
862 if (local_cov.size() > 0) {
863 for (
size_t i = 0,
nRows = local_cov.rows(),
nCols = local_cov.cols();
i <
nRows;
i++) {
864 for (
size_t j = 0; j <
nCols; ++j) {
865 v_local_cov.push_back(local_cov(
i, j));
871 (*m_CLhardware).push_back(
"STRIP");
881 (*m_CLparticleLink_eventIndex).push_back(particleLink_eventIndex);
882 (*m_CLparticleLink_barcode).push_back(particleLink_barcode);
883 (*m_CLbarcodesLinked).push_back(barcodesLinked);
884 (*m_CLparticle_charge).push_back(
charge);
885 (*m_CLetas).push_back(strip_ids);
886 (*m_CLphis).push_back(cst);
887 (*m_CLtots).push_back(tots);
905 (*m_CLlocal_cov).push_back(v_local_cov);
932 if (not xAODPixelSpacePointContainerHandle.isValid()) {
934 return StatusCode::FAILURE;
937 xAODPixelSPContainer = xAODPixelSpacePointContainerHandle.cptr();
942 if (not xAODStripSpacePointContainerHandle.isValid()) {
944 return StatusCode::FAILURE;
946 xAODStripSPContainer = xAODStripSpacePointContainerHandle.cptr();
951 if (not xAODStripSpacePointOverlapContainerHandle.isValid()) {
953 return StatusCode::FAILURE;
955 xAODStripSPOverlapContainer = xAODStripSpacePointOverlapContainerHandle.cptr();
960 if (xAODPixelSPContainer && xAODPixelSPContainer->
size() > 0) {
961 for (
const auto sp : *xAODPixelSPContainer) {
963 if (not linkAcc.isAvailable(*sp))
964 ATH_MSG_FATAL(
"no pixel SpacePoint link for xAOD::SpacePoint");
967 auto trk_sp = *linkAcc(*sp);
992 if (xAODStripSPContainer && xAODStripSPContainer->
size() > 0) {
995 for (
const auto sp : *xAODStripSPContainer) {
997 ATH_CHECK(striplinkAcc.isAvailable(*sp));
999 auto trk_sp = *striplinkAcc(*sp);
1019 std::vector<float> topstripDir(sp->topStripDirection().data(),
1020 sp->topStripDirection().data() +
1021 sp->topStripDirection().size());
1023 std::vector<float> botstripDir(sp->bottomStripDirection().data(),
1024 sp->bottomStripDirection().data() +
1025 sp->bottomStripDirection().size());
1027 std::vector<float> DstripCnt(sp->stripCenterDistance().data(),
1028 sp->stripCenterDistance().data() +
1029 sp->stripCenterDistance().size());
1031 std::vector<float> topstripCnt(sp->topStripCenter().data(),
1032 sp->topStripCenter().data() +
1033 sp->topStripCenter().size());
1035 (*m_SPtopStripDirection).push_back(topstripDir);
1036 (*m_SPbottomStripDirection).push_back(botstripDir);
1037 (*m_SPstripCenterDistance).push_back(DstripCnt);
1038 (*m_SPtopStripCenterPosition).push_back(topstripCnt);
1053 if (xAODStripSPOverlapContainer && xAODStripSPOverlapContainer->
size() > 0) {
1056 for (
const auto sp : *xAODStripSPOverlapContainer) {
1058 ATH_CHECK(stripOverlaplinkAcc.isAvailable(*sp));
1060 auto trk_sp = *stripOverlaplinkAcc(*sp);
1081 if ( flag<1 || flag > 3 )
1090 std::vector<float> topstripDir(sp->topStripDirection().data(),
1091 sp->topStripDirection().data() +
1092 sp->topStripDirection().size());
1094 std::vector<float> botstripDir(sp->bottomStripDirection().data(),
1095 sp->bottomStripDirection().data() +
1096 sp->bottomStripDirection().size());
1098 std::vector<float> DstripCnt(sp->stripCenterDistance().data(),
1099 sp->stripCenterDistance().data() +
1100 sp->stripCenterDistance().size());
1102 std::vector<float> topstripCnt(sp->topStripCenter().data(),
1103 sp->topStripCenter().data() +
1104 sp->topStripCenter().size());
1106 (*m_SPtopStripDirection).push_back(topstripDir);
1107 (*m_SPbottomStripDirection).push_back(botstripDir);
1108 (*m_SPstripCenterDistance).push_back(DstripCnt);
1109 (*m_SPtopStripCenterPosition).push_back(topstripCnt);
1129 if (not trackCollectionHandle.isValid()) {
1131 return StatusCode::FAILURE;
1133 trackCollection = trackCollectionHandle.cptr();
1137 if (not trackTruthCollectionHandle.isValid()) {
1139 return StatusCode::FAILURE;
1141 trackTruthCollection = trackTruthCollectionHandle.cptr();
1149 (*m_TRKproperties).clear();
1150 (*m_TRKpattern).clear();
1151 (*m_TRKperigee_position).clear();
1152 (*m_TRKperigee_momentum).clear();
1153 (*m_TRKmeasurementsOnTrack_pixcl_sctcl_index).clear();
1154 (*m_TRKoutliersOnTrack_pixcl_sctcl_index).clear();
1157 for (; trackIterator < (*trackCollection).end(); ++trackIterator) {
1158 if (!((*trackIterator))) {
1172 TrackTruthCollection::const_iterator
found = trackTruthCollection->find(tracklink2);
1174 const std::bitset<Trk::TrackInfo::NumberOfTrackProperties> &
properties =
info.properties();
1175 std::vector<int> v_properties;
1178 v_properties.push_back(
i);
1182 const std::bitset<Trk::TrackInfo::NumberOfTrackRecoInfo> &
pattern =
info.patternRecognition();
1183 std::vector<int> v_pattern;
1184 for (std::size_t
i = 0;
i <
pattern.size();
i++) {
1186 v_pattern.push_back(
i);
1196 std::vector<double> position,
momentum;
1207 position.push_back(0);
1208 position.push_back(0);
1209 position.push_back(0);
1216 if (measurementsOnTrack)
1217 mot = measurementsOnTrack->
size();
1218 if (outliersOnTrack)
1219 oot = outliersOnTrack->
size();
1220 std::vector<int> measurementsOnTrack_pixcl_sctcl_index, outliersOnTrack_pixcl_sctcl_index;
1221 int TTCindex, TTCevent_index, TTCparticle_link;
1222 float TTCprobability;
1223 if (measurementsOnTrack) {
1224 for (
size_t i = 0;
i < measurementsOnTrack->
size();
i++) {
1229 measurementsOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[pixcl->
prepRawData()->
identify()]);
1232 measurementsOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[sctcl->
prepRawData()->
identify()]);
1234 measurementsOnTrack_pixcl_sctcl_index.push_back(-1);
1238 if (outliersOnTrack) {
1239 for (
size_t i = 0;
i < outliersOnTrack->
size();
i++) {
1244 outliersOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[pixcl->
prepRawData()->
identify()]);
1246 outliersOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[sctcl->
prepRawData()->
identify()]);
1248 outliersOnTrack_pixcl_sctcl_index.push_back(-1);
1252 if (
found != trackTruthCollection->end()) {
1253 TTCindex =
found->first.index();
1254 TTCevent_index =
found->second.particleLink().eventIndex();
1255 TTCparticle_link =
found->second.particleLink().barcode();
1256 TTCprobability =
found->second.probability();
1258 TTCindex = TTCevent_index = TTCparticle_link = -999;
1259 TTCprobability = -1;
1267 (*m_TRKproperties).push_back(v_properties);
1268 (*m_TRKpattern).push_back(v_pattern);
1271 (*m_TRKmeasurementsOnTrack_pixcl_sctcl_index).push_back(measurementsOnTrack_pixcl_sctcl_index);
1272 (*m_TRKoutliersOnTrack_pixcl_sctcl_index).push_back(outliersOnTrack_pixcl_sctcl_index);
1274 (*m_TRKperigee_position).push_back(position);
1275 (*m_TRKperigee_momentum).push_back(
momentum);
1295 if (not detailedTrackTruthCollectionHandle.isValid()) {
1297 return StatusCode::FAILURE;
1299 detailedTrackTruthCollection = detailedTrackTruthCollectionHandle.cptr();
1303 (*m_DTTtrajectory_eventindex).clear();
1304 (*m_DTTtrajectory_barcode).clear();
1305 (*m_DTTstTruth_subDetType).clear();
1306 (*m_DTTstTrack_subDetType).clear();
1307 (*m_DTTstCommon_subDetType).clear();
1311 DetailedTrackTruthCollection::const_iterator detailedTrackTruthIterator = (*detailedTrackTruthCollection).begin();
1312 for (; detailedTrackTruthIterator != (*detailedTrackTruthCollection).end(); ++detailedTrackTruthIterator) {
1313 std::vector<int> DTTtrajectory_eventindex, DTTtrajectory_barcode, DTTstTruth_subDetType, DTTstTrack_subDetType,
1314 DTTstCommon_subDetType;
1315 const TruthTrajectory &traj = detailedTrackTruthIterator->second.trajectory();
1316 for (
size_t j = 0; j < traj.size(); j++) {
1317 DTTtrajectory_eventindex.push_back(traj[j].
eventIndex());
1318 DTTtrajectory_barcode.push_back(traj[j].
barcode());
1336 (*m_DTTtrajectory_eventindex).push_back(DTTtrajectory_eventindex);
1337 (*m_DTTtrajectory_barcode).push_back(DTTtrajectory_barcode);
1338 (*m_DTTstTruth_subDetType).push_back(DTTstTruth_subDetType);
1339 (*m_DTTstTrack_subDetType).push_back(DTTstTrack_subDetType);
1340 (*m_DTTstCommon_subDetType).push_back(DTTstCommon_subDetType);
1351 return StatusCode::SUCCESS;
1.8.18