376 const EventContext &ctx = Gaudi::Hive::currentContext();
382 std::map<Identifier, long int> clusterIDMapIdx;
385 std::map<Identifier, long int> clusterIDMapSpacePointIdx;
389 std::map<std::pair<int, int>, std::pair<bool, int>> allTruthParticles;
393 if (not mcEventCollectionHandle.isValid()) {
395 return StatusCode::FAILURE;
397 mcCollptr = mcEventCollectionHandle.cptr();
402 if (not eventInfoHandle.isValid()) {
404 return StatusCode::FAILURE;
406 eventInfo = eventInfoHandle.cptr();
411 std::map<int, int> allSubEvents;
415 bool duplicateSubeventID =
false;
416 for (
unsigned int cntr = 0; cntr < mcCollptr->
size(); ++cntr) {
417 int ID = mcCollptr->
at(cntr)->event_number();
426 if (
it == allSubEvents.end())
427 allSubEvents.insert(std::make_pair(
ID, 1));
430 duplicateSubeventID =
true;
434 if (duplicateSubeventID) {
441 (*m_Part_vParentID).clear();
442 (*m_Part_vParentBarcode).clear();
445 for (
unsigned int cntr = 0; cntr < mcCollptr->
size(); ++cntr) {
446 const HepMC::GenEvent *genEvt = (mcCollptr->
at(cntr));
455 for (
auto p : *genEvt) {
457 float px,
py,
pz,
pt,
eta, vx, vy, vz,
radius,
status,
charge = 0.;
458 std::vector<int> vParentID;
459 std::vector<int> vParentBarcode;
461 int vProdNin, vProdNout, vProdStatus, vProdBarcode;
462 bool passed =
isPassed(
p,
px,
py,
pz,
pt, eta, vx, vy, vz,
radius,
status,
charge, vParentID, vParentBarcode,
463 vProdNin, vProdNout, vProdStatus, vProdBarcode);
464 allTruthParticles.insert(std::make_pair(std::make_pair(genEvt->event_number(),
HepMC::barcode(
p)),
465 std::make_pair(
passed, 0)));
487 (*m_Part_vParentID).push_back(vParentID);
488 (*m_Part_vParentBarcode).push_back(vParentBarcode);
501 if (not pixelClusterContainerHandle.isValid()) {
503 return StatusCode::FAILURE;
509 if (not stripClusterContainerHandle.isValid()) {
511 return StatusCode::FAILURE;
515 auto cartesion_to_spherical = [](
const Amg::Vector3D &xyzVec,
float &eta_,
float &phi_) {
518 r3 += xyzVec[
idx] * xyzVec[
idx];
521 phi_ = atan2(xyzVec[1], xyzVec[0]);
522 float theta_ = acos(xyzVec[2] / r3);
523 eta_ =
log(
tan(0.5 * theta_));
532 (*m_CLhardware).clear();
533 (*m_CLparticleLink_eventIndex).clear();
534 (*m_CLparticleLink_barcode).clear();
535 (*m_CLbarcodesLinked).clear();
536 (*m_CLparticle_charge).clear();
540 (*m_CLlocal_cov).clear();
547 if (not sdoCollectionHandle.isValid()) {
549 return StatusCode::FAILURE;
551 sdoCollection = sdoCollectionHandle.cptr();
555 if (clusterCollection->empty())
566 float norm_x = fabs(my_normal.x()) > 1
e-5 ? my_normal.x() : 0.;
567 float norm_y = fabs(my_normal.y()) > 1
e-5 ? my_normal.y() : 0.;
568 float norm_z = fabs(my_normal.z()) > 1
e-5 ? my_normal.z() : 0.;
573 ATH_MSG_ERROR(
"Dynamic cast failed at " << __LINE__ <<
" of MergedPixelsTool.cxx.");
574 return StatusCode::FAILURE;
578 for (
const auto cluster : *clusterCollection) {
584 const Amg::MatrixX &local_cov = cluster->localCovariance();
586 std::vector<std::pair<int, int>>
barcodes = {};
587 std::vector<int> particleLink_eventIndex = {};
588 std::vector<int> particleLink_barcode = {};
589 std::vector<bool> barcodesLinked = {};
590 std::vector<float>
charge = {};
591 std::vector<int> phis = {};
592 std::vector<int>
etas = {};
593 std::vector<int> tots = {};
599 float charge_count = 0;
602 for (
unsigned int rdo = 0; rdo < cluster->rdoList().
size(); rdo++) {
603 const auto &rdoID = cluster->rdoList().at(rdo);
616 charge_count += cluster->totList().at(rdo);
620 tots.push_back(cluster->totList().at(rdo));
622 auto pos = sdoCollection->find(rdoID);
623 if (
pos != sdoCollection->end()) {
624 for (
auto deposit :
pos->second.getdeposits()) {
630 particleLink_eventIndex.push_back(particleLink.
eventIndex());
631 particleLink_barcode.push_back(particleLink.
barcode());
632 charge.push_back(deposit.second);
633 barcodesLinked.push_back(particleLink.
isValid());
650 Amg::Vector3D localDirection = localEndPosition - localStartPosition;
652 float loc_eta = 0, loc_phi = 0;
653 cartesion_to_spherical(localDirection, loc_eta, loc_phi);
658 Amg::Vector3D direction = globalEndPosition - globalStartPosition;
659 float glob_eta = 0, glob_phi = 0;
660 cartesion_to_spherical(direction, glob_eta, glob_phi);
665 float trkphicomp = direction.dot(my_phiax);
666 float trketacomp = direction.dot(my_etaax);
667 float trknormcomp = direction.dot(my_normal);
668 double phi_angle = atan2(trknormcomp, trkphicomp);
669 double eta_angle = atan2(trknormcomp, trketacomp);
671 clusterIDMapIdx[cluster->identify()] =
m_selected;
672 std::vector<double> v_local_cov;
673 if (local_cov.size() > 0) {
674 for (
size_t i = 0,
nRows = local_cov.rows(),
nCols = local_cov.cols();
i <
nRows;
i++) {
675 for (
size_t j = 0; j <
nCols; ++j) {
676 v_local_cov.push_back(local_cov(
i, j));
683 (*m_CLhardware).push_back(
"PIXEL");
693 (*m_CLparticleLink_eventIndex).push_back(particleLink_eventIndex);
694 (*m_CLparticleLink_barcode).push_back(particleLink_barcode);
695 (*m_CLbarcodesLinked).push_back(barcodesLinked);
696 (*m_CLparticle_charge).push_back(
charge);
697 (*m_CLetas).push_back(
etas);
698 (*m_CLphis).push_back(phis);
699 (*m_CLtots).push_back(tots);
717 (*m_CLlocal_cov).push_back(v_local_cov);
736 if (not sdoCollectionHandle.isValid()) {
738 return StatusCode::FAILURE;
740 sdoCollection = sdoCollectionHandle.cptr();
744 if (clusterCollection->empty())
756 float norm_x = fabs(my_normal.x()) > 1
e-5 ? my_normal.x() : 0.;
757 float norm_y = fabs(my_normal.y()) > 1
e-5 ? my_normal.y() : 0.;
758 float norm_z = fabs(my_normal.z()) > 1
e-5 ? my_normal.z() : 0.;
761 for (
const auto cluster : *clusterCollection) {
767 const Amg::MatrixX &local_cov = cluster->localCovariance();
769 std::vector<std::pair<int, int>>
barcodes = {};
770 std::vector<int> particleLink_eventIndex = {};
771 std::vector<int> particleLink_barcode = {};
772 std::vector<bool> barcodesLinked = {};
773 std::vector<float>
charge = {};
775 std::vector<int> tots = {};
776 std::vector<int> strip_ids = {};
778 int max_strip = -999;
780 float charge_count = 0;
783 for (
unsigned int rdo = 0; rdo < cluster->rdoList().
size(); rdo++) {
784 const auto &rdoID = cluster->rdoList().at(rdo);
788 if (min_strip > strip)
790 if (max_strip < strip)
792 strip_ids.push_back(strip);
797 auto pos = sdoCollection->find(rdoID);
798 if (
pos != sdoCollection->end()) {
799 for (
auto deposit :
pos->second.getdeposits()) {
806 particleLink_eventIndex.push_back(particleLink.
eventIndex());
807 particleLink_barcode.push_back(particleLink.
barcode());
808 charge.push_back(deposit.second);
809 barcodesLinked.push_back(particleLink.
isValid());
819 ATH_MSG_ERROR(
"Failed at " << __LINE__ <<
" of accessing SCT ModuleSide Design");
820 return StatusCode::FAILURE;
824 std::pair<Amg::Vector3D, Amg::Vector3D> ends(
838 Amg::Vector3D localDirection = localEndPosition - localStartPosition;
839 float loc_eta = 0, loc_phi = 0;
840 cartesion_to_spherical(localDirection, loc_eta, loc_phi);
845 Amg::Vector3D direction = globalEndPosition - globalStartPosition;
846 float glob_eta = 0, glob_phi = 0;
847 cartesion_to_spherical(direction, glob_eta, glob_phi);
852 float trkphicomp = direction.dot(my_phiax);
853 float trketacomp = direction.dot(my_etaax);
854 float trknormcomp = direction.dot(my_normal);
855 double phi_angle = atan2(trknormcomp, trkphicomp);
856 double eta_angle = atan2(trknormcomp, trketacomp);
859 clusterIDMapIdx[cluster->identify()] =
m_selected;
861 std::vector<int> cst;
862 for (
unsigned strip = 0;
strip < strip_ids.size();
strip++) {
865 std::vector<double> v_local_cov;
866 if (local_cov.size() > 0) {
867 for (
size_t i = 0,
nRows = local_cov.rows(),
nCols = local_cov.cols();
i <
nRows;
i++) {
868 for (
size_t j = 0; j <
nCols; ++j) {
869 v_local_cov.push_back(local_cov(
i, j));
875 (*m_CLhardware).push_back(
"STRIP");
885 (*m_CLparticleLink_eventIndex).push_back(particleLink_eventIndex);
886 (*m_CLparticleLink_barcode).push_back(particleLink_barcode);
887 (*m_CLbarcodesLinked).push_back(barcodesLinked);
888 (*m_CLparticle_charge).push_back(
charge);
889 (*m_CLetas).push_back(strip_ids);
890 (*m_CLphis).push_back(cst);
891 (*m_CLtots).push_back(tots);
909 (*m_CLlocal_cov).push_back(v_local_cov);
936 if (not xAODPixelSpacePointContainerHandle.isValid()) {
938 return StatusCode::FAILURE;
941 xAODPixelSPContainer = xAODPixelSpacePointContainerHandle.cptr();
946 if (not xAODStripSpacePointContainerHandle.isValid()) {
948 return StatusCode::FAILURE;
950 xAODStripSPContainer = xAODStripSpacePointContainerHandle.cptr();
955 if (not xAODStripSpacePointOverlapContainerHandle.isValid()) {
957 return StatusCode::FAILURE;
959 xAODStripSPOverlapContainer = xAODStripSpacePointOverlapContainerHandle.cptr();
964 if (xAODPixelSPContainer && xAODPixelSPContainer->
size() > 0) {
965 for (
const auto sp : *xAODPixelSPContainer) {
967 if (not linkAcc.isAvailable(*sp))
968 ATH_MSG_FATAL(
"no pixel SpacePoint link for xAOD::SpacePoint");
971 auto trk_sp = *linkAcc(*sp);
996 if (xAODStripSPContainer && xAODStripSPContainer->
size() > 0) {
999 for (
const auto sp : *xAODStripSPContainer) {
1001 ATH_CHECK(striplinkAcc.isAvailable(*sp));
1003 auto trk_sp = *striplinkAcc(*sp);
1023 std::vector<float> topstripDir(sp->topStripDirection().data(),
1024 sp->topStripDirection().data() +
1025 sp->topStripDirection().size());
1027 std::vector<float> botstripDir(sp->bottomStripDirection().data(),
1028 sp->bottomStripDirection().data() +
1029 sp->bottomStripDirection().size());
1031 std::vector<float> DstripCnt(sp->stripCenterDistance().data(),
1032 sp->stripCenterDistance().data() +
1033 sp->stripCenterDistance().size());
1035 std::vector<float> topstripCnt(sp->topStripCenter().data(),
1036 sp->topStripCenter().data() +
1037 sp->topStripCenter().size());
1039 (*m_SPtopStripDirection).push_back(topstripDir);
1040 (*m_SPbottomStripDirection).push_back(botstripDir);
1041 (*m_SPstripCenterDistance).push_back(DstripCnt);
1042 (*m_SPtopStripCenterPosition).push_back(topstripCnt);
1057 if (xAODStripSPOverlapContainer && xAODStripSPOverlapContainer->
size() > 0) {
1060 for (
const auto sp : *xAODStripSPOverlapContainer) {
1062 ATH_CHECK(stripOverlaplinkAcc.isAvailable(*sp));
1064 auto trk_sp = *stripOverlaplinkAcc(*sp);
1085 if ( flag<1 || flag > 2 )
1094 std::vector<float> topstripDir(sp->topStripDirection().data(),
1095 sp->topStripDirection().data() +
1096 sp->topStripDirection().size());
1098 std::vector<float> botstripDir(sp->bottomStripDirection().data(),
1099 sp->bottomStripDirection().data() +
1100 sp->bottomStripDirection().size());
1102 std::vector<float> DstripCnt(sp->stripCenterDistance().data(),
1103 sp->stripCenterDistance().data() +
1104 sp->stripCenterDistance().size());
1106 std::vector<float> topstripCnt(sp->topStripCenter().data(),
1107 sp->topStripCenter().data() +
1108 sp->topStripCenter().size());
1110 (*m_SPtopStripDirection).push_back(topstripDir);
1111 (*m_SPbottomStripDirection).push_back(botstripDir);
1112 (*m_SPstripCenterDistance).push_back(DstripCnt);
1113 (*m_SPtopStripCenterPosition).push_back(topstripCnt);
1133 if (not trackCollectionHandle.isValid()) {
1135 return StatusCode::FAILURE;
1137 trackCollection = trackCollectionHandle.cptr();
1141 if (not trackTruthCollectionHandle.isValid()) {
1143 return StatusCode::FAILURE;
1145 trackTruthCollection = trackTruthCollectionHandle.cptr();
1153 (*m_TRKproperties).clear();
1154 (*m_TRKpattern).clear();
1155 (*m_TRKperigee_position).clear();
1156 (*m_TRKperigee_momentum).clear();
1157 (*m_TRKmeasurementsOnTrack_pixcl_sctcl_index).clear();
1158 (*m_TRKoutliersOnTrack_pixcl_sctcl_index).clear();
1161 for (; trackIterator < (*trackCollection).end(); ++trackIterator) {
1162 if (!((*trackIterator))) {
1176 TrackTruthCollection::const_iterator
found = trackTruthCollection->find(tracklink2);
1178 const std::bitset<Trk::TrackInfo::NumberOfTrackProperties> &
properties =
info.properties();
1179 std::vector<int> v_properties;
1182 v_properties.push_back(
i);
1186 const std::bitset<Trk::TrackInfo::NumberOfTrackRecoInfo> &
pattern =
info.patternRecognition();
1187 std::vector<int> v_pattern;
1188 for (std::size_t
i = 0;
i <
pattern.size();
i++) {
1190 v_pattern.push_back(
i);
1200 std::vector<double> position,
momentum;
1211 position.push_back(0);
1212 position.push_back(0);
1213 position.push_back(0);
1220 if (measurementsOnTrack)
1221 mot = measurementsOnTrack->
size();
1222 if (outliersOnTrack)
1223 oot = outliersOnTrack->
size();
1224 std::vector<int> measurementsOnTrack_pixcl_sctcl_index, outliersOnTrack_pixcl_sctcl_index;
1225 int TTCindex, TTCevent_index, TTCparticle_link;
1226 float TTCprobability;
1227 if (measurementsOnTrack) {
1228 for (
size_t i = 0;
i < measurementsOnTrack->
size();
i++) {
1233 measurementsOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[pixcl->
prepRawData()->
identify()]);
1236 measurementsOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[sctcl->
prepRawData()->
identify()]);
1238 measurementsOnTrack_pixcl_sctcl_index.push_back(-1);
1242 if (outliersOnTrack) {
1243 for (
size_t i = 0;
i < outliersOnTrack->
size();
i++) {
1248 outliersOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[pixcl->
prepRawData()->
identify()]);
1250 outliersOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[sctcl->
prepRawData()->
identify()]);
1252 outliersOnTrack_pixcl_sctcl_index.push_back(-1);
1256 if (
found != trackTruthCollection->end()) {
1257 TTCindex =
found->first.index();
1258 TTCevent_index =
found->second.particleLink().eventIndex();
1259 TTCparticle_link =
found->second.particleLink().barcode();
1260 TTCprobability =
found->second.probability();
1262 TTCindex = TTCevent_index = TTCparticle_link = -999;
1263 TTCprobability = -1;
1271 (*m_TRKproperties).push_back(v_properties);
1272 (*m_TRKpattern).push_back(v_pattern);
1275 (*m_TRKmeasurementsOnTrack_pixcl_sctcl_index).push_back(measurementsOnTrack_pixcl_sctcl_index);
1276 (*m_TRKoutliersOnTrack_pixcl_sctcl_index).push_back(outliersOnTrack_pixcl_sctcl_index);
1278 (*m_TRKperigee_position).push_back(position);
1279 (*m_TRKperigee_momentum).push_back(
momentum);
1299 if (not detailedTrackTruthCollectionHandle.isValid()) {
1301 return StatusCode::FAILURE;
1303 detailedTrackTruthCollection = detailedTrackTruthCollectionHandle.cptr();
1307 (*m_DTTtrajectory_eventindex).clear();
1308 (*m_DTTtrajectory_barcode).clear();
1309 (*m_DTTstTruth_subDetType).clear();
1310 (*m_DTTstTrack_subDetType).clear();
1311 (*m_DTTstCommon_subDetType).clear();
1315 DetailedTrackTruthCollection::const_iterator detailedTrackTruthIterator = (*detailedTrackTruthCollection).begin();
1316 for (; detailedTrackTruthIterator != (*detailedTrackTruthCollection).end(); ++detailedTrackTruthIterator) {
1317 std::vector<int> DTTtrajectory_eventindex, DTTtrajectory_barcode, DTTstTruth_subDetType, DTTstTrack_subDetType,
1318 DTTstCommon_subDetType;
1319 const TruthTrajectory &traj = detailedTrackTruthIterator->second.trajectory();
1320 for (
size_t j = 0; j < traj.size(); j++) {
1321 DTTtrajectory_eventindex.push_back(traj[j].
eventIndex());
1322 DTTtrajectory_barcode.push_back(traj[j].
barcode());
1340 (*m_DTTtrajectory_eventindex).push_back(DTTtrajectory_eventindex);
1341 (*m_DTTtrajectory_barcode).push_back(DTTtrajectory_barcode);
1342 (*m_DTTstTruth_subDetType).push_back(DTTstTruth_subDetType);
1343 (*m_DTTstTrack_subDetType).push_back(DTTstTrack_subDetType);
1344 (*m_DTTstCommon_subDetType).push_back(DTTstCommon_subDetType);
1355 return StatusCode::SUCCESS;
1.8.18