375 const EventContext &ctx = Gaudi::Hive::currentContext();
381 std::map<Identifier, long int> clusterIDMapIdx;
384 std::map<Identifier, long int> clusterIDMapSpacePointIdx;
388 std::map<std::pair<int, int>, std::pair<bool, int>> allTruthParticles;
392 if (not mcEventCollectionHandle.isValid()) {
394 return StatusCode::FAILURE;
396 mcCollptr = mcEventCollectionHandle.cptr();
401 if (not eventInfoHandle.isValid()) {
403 return StatusCode::FAILURE;
410 std::map<int, int> allSubEvents;
414 bool duplicateSubeventID =
false;
415 for (
unsigned int cntr = 0; cntr < mcCollptr->
size(); ++cntr) {
416 int ID = mcCollptr->
at(cntr)->event_number();
425 if (
it == allSubEvents.end())
426 allSubEvents.insert(std::make_pair(
ID, 1));
429 duplicateSubeventID =
true;
433 if (duplicateSubeventID) {
440 (*m_Part_vParentID).clear();
441 (*m_Part_vParentBarcode).clear();
444 for (
unsigned int cntr = 0; cntr < mcCollptr->
size(); ++cntr) {
445 const HepMC::GenEvent *genEvt = (mcCollptr->
at(cntr));
454 for (
auto p : *genEvt) {
456 float px,
py,
pz,
pt,
eta, vx, vy, vz,
radius,
status,
charge = 0.;
457 std::vector<int> vParentID;
458 std::vector<int> vParentBarcode;
460 int vProdNin, vProdNout, vProdStatus, vProdBarcode;
461 bool passed =
isPassed(
p,
px,
py,
pz,
pt, eta, vx, vy, vz,
radius,
status,
charge, vParentID, vParentBarcode,
462 vProdNin, vProdNout, vProdStatus, vProdBarcode);
463 allTruthParticles.insert(std::make_pair(std::make_pair(genEvt->event_number(),
HepMC::barcode(
p)),
464 std::make_pair(
passed, 0)));
486 (*m_Part_vParentID).push_back(vParentID);
487 (*m_Part_vParentBarcode).push_back(vParentBarcode);
500 if (not pixelClusterContainerHandle.isValid()) {
502 return StatusCode::FAILURE;
508 if (not stripClusterContainerHandle.isValid()) {
510 return StatusCode::FAILURE;
514 auto cartesion_to_spherical = [](
const Amg::Vector3D &xyzVec,
float &eta_,
float &phi_) {
517 r3 += xyzVec[
idx] * xyzVec[
idx];
520 phi_ = atan2(xyzVec[1], xyzVec[0]);
521 float theta_ = acos(xyzVec[2] / r3);
522 eta_ =
log(
tan(0.5 * theta_));
531 (*m_CLhardware).clear();
532 (*m_CLparticleLink_eventIndex).clear();
533 (*m_CLparticleLink_barcode).clear();
534 (*m_CLbarcodesLinked).clear();
535 (*m_CLparticle_charge).clear();
539 (*m_CLlocal_cov).clear();
546 if (not sdoCollectionHandle.isValid()) {
548 return StatusCode::FAILURE;
550 sdoCollection = sdoCollectionHandle.cptr();
554 if (clusterCollection->empty())
565 float norm_x = fabs(my_normal.x()) > 1
e-5 ? my_normal.x() : 0.;
566 float norm_y = fabs(my_normal.y()) > 1
e-5 ? my_normal.y() : 0.;
567 float norm_z = fabs(my_normal.z()) > 1
e-5 ? my_normal.z() : 0.;
572 ATH_MSG_ERROR(
"Dynamic cast failed at " << __LINE__ <<
" of MergedPixelsTool.cxx.");
573 return StatusCode::FAILURE;
577 for (
const auto cluster : *clusterCollection) {
585 std::vector<std::pair<int, int>>
barcodes = {};
586 std::vector<int> particleLink_eventIndex = {};
587 std::vector<int> particleLink_barcode = {};
588 std::vector<bool> barcodesLinked = {};
589 std::vector<float>
charge = {};
590 std::vector<int> phis = {};
591 std::vector<int>
etas = {};
592 std::vector<int> tots = {};
598 float charge_count = 0;
601 for (
unsigned int rdo = 0; rdo <
cluster->rdoList().
size(); rdo++) {
602 const auto &rdoID =
cluster->rdoList().at(rdo);
615 charge_count +=
cluster->totList().at(rdo);
619 tots.push_back(
cluster->totList().at(rdo));
621 auto pos = sdoCollection->find(rdoID);
622 if (
pos != sdoCollection->end()) {
623 for (
auto deposit :
pos->second.getdeposits()) {
629 particleLink_eventIndex.push_back(particleLink.
eventIndex());
630 particleLink_barcode.push_back(particleLink.
barcode());
631 charge.push_back(deposit.second);
632 barcodesLinked.push_back(particleLink.
isValid());
649 Amg::Vector3D localDirection = localEndPosition - localStartPosition;
651 float loc_eta = 0, loc_phi = 0;
652 cartesion_to_spherical(localDirection, loc_eta, loc_phi);
657 Amg::Vector3D direction = globalEndPosition - globalStartPosition;
658 float glob_eta = 0, glob_phi = 0;
659 cartesion_to_spherical(direction, glob_eta, glob_phi);
664 float trkphicomp = direction.dot(my_phiax);
665 float trketacomp = direction.dot(my_etaax);
666 float trknormcomp = direction.dot(my_normal);
667 double phi_angle = atan2(trknormcomp, trkphicomp);
668 double eta_angle = atan2(trknormcomp, trketacomp);
671 std::vector<double> v_local_cov;
672 if (local_cov.size() > 0) {
673 for (
size_t i = 0,
nRows = local_cov.rows(),
nCols = local_cov.cols();
i <
nRows;
i++) {
674 for (
size_t j = 0; j <
nCols; ++j) {
675 v_local_cov.push_back(local_cov(
i, j));
682 (*m_CLhardware).push_back(
"PIXEL");
692 (*m_CLparticleLink_eventIndex).push_back(particleLink_eventIndex);
693 (*m_CLparticleLink_barcode).push_back(particleLink_barcode);
694 (*m_CLbarcodesLinked).push_back(barcodesLinked);
695 (*m_CLparticle_charge).push_back(
charge);
696 (*m_CLetas).push_back(
etas);
697 (*m_CLphis).push_back(phis);
698 (*m_CLtots).push_back(tots);
716 (*m_CLlocal_cov).push_back(v_local_cov);
735 if (not sdoCollectionHandle.isValid()) {
737 return StatusCode::FAILURE;
739 sdoCollection = sdoCollectionHandle.cptr();
743 if (clusterCollection->empty())
755 float norm_x = fabs(my_normal.x()) > 1
e-5 ? my_normal.x() : 0.;
756 float norm_y = fabs(my_normal.y()) > 1
e-5 ? my_normal.y() : 0.;
757 float norm_z = fabs(my_normal.z()) > 1
e-5 ? my_normal.z() : 0.;
760 for (
const auto cluster : *clusterCollection) {
768 std::vector<std::pair<int, int>>
barcodes = {};
769 std::vector<int> particleLink_eventIndex = {};
770 std::vector<int> particleLink_barcode = {};
771 std::vector<bool> barcodesLinked = {};
772 std::vector<float>
charge = {};
774 std::vector<int> tots = {};
775 std::vector<int> strip_ids = {};
777 int max_strip = -999;
779 float charge_count = 0;
782 for (
unsigned int rdo = 0; rdo <
cluster->rdoList().
size(); rdo++) {
783 const auto &rdoID =
cluster->rdoList().at(rdo);
787 if (min_strip > strip)
789 if (max_strip < strip)
791 strip_ids.push_back(strip);
796 auto pos = sdoCollection->find(rdoID);
797 if (
pos != sdoCollection->end()) {
798 for (
auto deposit :
pos->second.getdeposits()) {
805 particleLink_eventIndex.push_back(particleLink.
eventIndex());
806 particleLink_barcode.push_back(particleLink.
barcode());
807 charge.push_back(deposit.second);
808 barcodesLinked.push_back(particleLink.
isValid());
818 ATH_MSG_ERROR(
"Failed at " << __LINE__ <<
" of accessing SCT ModuleSide Design");
819 return StatusCode::FAILURE;
823 std::pair<Amg::Vector3D, Amg::Vector3D> ends(
837 Amg::Vector3D localDirection = localEndPosition - localStartPosition;
838 float loc_eta = 0, loc_phi = 0;
839 cartesion_to_spherical(localDirection, loc_eta, loc_phi);
844 Amg::Vector3D direction = globalEndPosition - globalStartPosition;
845 float glob_eta = 0, glob_phi = 0;
846 cartesion_to_spherical(direction, glob_eta, glob_phi);
851 float trkphicomp = direction.dot(my_phiax);
852 float trketacomp = direction.dot(my_etaax);
853 float trknormcomp = direction.dot(my_normal);
854 double phi_angle = atan2(trknormcomp, trkphicomp);
855 double eta_angle = atan2(trknormcomp, trketacomp);
860 std::vector<int> cst;
861 for (
unsigned strip = 0;
strip < strip_ids.size();
strip++) {
864 std::vector<double> v_local_cov;
865 if (local_cov.size() > 0) {
866 for (
size_t i = 0,
nRows = local_cov.rows(),
nCols = local_cov.cols();
i <
nRows;
i++) {
867 for (
size_t j = 0; j <
nCols; ++j) {
868 v_local_cov.push_back(local_cov(
i, j));
874 (*m_CLhardware).push_back(
"STRIP");
884 (*m_CLparticleLink_eventIndex).push_back(particleLink_eventIndex);
885 (*m_CLparticleLink_barcode).push_back(particleLink_barcode);
886 (*m_CLbarcodesLinked).push_back(barcodesLinked);
887 (*m_CLparticle_charge).push_back(
charge);
888 (*m_CLetas).push_back(strip_ids);
889 (*m_CLphis).push_back(cst);
890 (*m_CLtots).push_back(tots);
908 (*m_CLlocal_cov).push_back(v_local_cov);
935 if (not xAODPixelSpacePointContainerHandle.isValid()) {
937 return StatusCode::FAILURE;
940 xAODPixelSPContainer = xAODPixelSpacePointContainerHandle.cptr();
945 if (not xAODStripSpacePointContainerHandle.isValid()) {
947 return StatusCode::FAILURE;
949 xAODStripSPContainer = xAODStripSpacePointContainerHandle.cptr();
954 if (not xAODStripSpacePointOverlapContainerHandle.isValid()) {
956 return StatusCode::FAILURE;
958 xAODStripSPOverlapContainer = xAODStripSpacePointOverlapContainerHandle.cptr();
963 if (xAODPixelSPContainer && xAODPixelSPContainer->
size() > 0) {
964 for (
const auto sp : *xAODPixelSPContainer) {
966 if (not linkAcc.isAvailable(*sp))
967 ATH_MSG_FATAL(
"no pixel SpacePoint link for xAOD::SpacePoint");
970 auto trk_sp = *linkAcc(*sp);
995 if (xAODStripSPContainer && xAODStripSPContainer->
size() > 0) {
998 for (
const auto sp : *xAODStripSPContainer) {
1000 ATH_CHECK(striplinkAcc.isAvailable(*sp));
1002 auto trk_sp = *striplinkAcc(*sp);
1022 std::vector<float> topstripDir(sp->topStripDirection().data(),
1023 sp->topStripDirection().data() +
1024 sp->topStripDirection().size());
1026 std::vector<float> botstripDir(sp->bottomStripDirection().data(),
1027 sp->bottomStripDirection().data() +
1028 sp->bottomStripDirection().size());
1030 std::vector<float> DstripCnt(sp->stripCenterDistance().data(),
1031 sp->stripCenterDistance().data() +
1032 sp->stripCenterDistance().size());
1034 std::vector<float> topstripCnt(sp->topStripCenter().data(),
1035 sp->topStripCenter().data() +
1036 sp->topStripCenter().size());
1038 (*m_SPtopStripDirection).push_back(topstripDir);
1039 (*m_SPbottomStripDirection).push_back(botstripDir);
1040 (*m_SPstripCenterDistance).push_back(DstripCnt);
1041 (*m_SPtopStripCenterPosition).push_back(topstripCnt);
1056 if (xAODStripSPOverlapContainer && xAODStripSPOverlapContainer->
size() > 0) {
1059 for (
const auto sp : *xAODStripSPOverlapContainer) {
1061 ATH_CHECK(stripOverlaplinkAcc.isAvailable(*sp));
1063 auto trk_sp = *stripOverlaplinkAcc(*sp);
1084 if ( flag<1 || flag > 3 )
1093 std::vector<float> topstripDir(sp->topStripDirection().data(),
1094 sp->topStripDirection().data() +
1095 sp->topStripDirection().size());
1097 std::vector<float> botstripDir(sp->bottomStripDirection().data(),
1098 sp->bottomStripDirection().data() +
1099 sp->bottomStripDirection().size());
1101 std::vector<float> DstripCnt(sp->stripCenterDistance().data(),
1102 sp->stripCenterDistance().data() +
1103 sp->stripCenterDistance().size());
1105 std::vector<float> topstripCnt(sp->topStripCenter().data(),
1106 sp->topStripCenter().data() +
1107 sp->topStripCenter().size());
1109 (*m_SPtopStripDirection).push_back(topstripDir);
1110 (*m_SPbottomStripDirection).push_back(botstripDir);
1111 (*m_SPstripCenterDistance).push_back(DstripCnt);
1112 (*m_SPtopStripCenterPosition).push_back(topstripCnt);
1132 if (not trackCollectionHandle.isValid()) {
1134 return StatusCode::FAILURE;
1136 trackCollection = trackCollectionHandle.cptr();
1140 if (not trackTruthCollectionHandle.isValid()) {
1142 return StatusCode::FAILURE;
1144 trackTruthCollection = trackTruthCollectionHandle.cptr();
1152 (*m_TRKproperties).clear();
1153 (*m_TRKpattern).clear();
1154 (*m_TRKperigee_position).clear();
1155 (*m_TRKperigee_momentum).clear();
1156 (*m_TRKmeasurementsOnTrack_pixcl_sctcl_index).clear();
1157 (*m_TRKoutliersOnTrack_pixcl_sctcl_index).clear();
1160 for (; trackIterator < (*trackCollection).end(); ++trackIterator) {
1161 if (!((*trackIterator))) {
1175 TrackTruthCollection::const_iterator
found = trackTruthCollection->find(tracklink2);
1177 const std::bitset<Trk::TrackInfo::NumberOfTrackProperties> &
properties =
info.properties();
1178 std::vector<int> v_properties;
1181 v_properties.push_back(
i);
1185 const std::bitset<Trk::TrackInfo::NumberOfTrackRecoInfo> &
pattern =
info.patternRecognition();
1186 std::vector<int> v_pattern;
1187 for (std::size_t
i = 0;
i <
pattern.size();
i++) {
1189 v_pattern.push_back(
i);
1199 std::vector<double> position,
momentum;
1210 position.push_back(0);
1211 position.push_back(0);
1212 position.push_back(0);
1219 if (measurementsOnTrack)
1220 mot = measurementsOnTrack->
size();
1221 if (outliersOnTrack)
1222 oot = outliersOnTrack->
size();
1223 std::vector<int> measurementsOnTrack_pixcl_sctcl_index, outliersOnTrack_pixcl_sctcl_index;
1224 int TTCindex, TTCevent_index, TTCparticle_link;
1225 float TTCprobability;
1226 if (measurementsOnTrack) {
1227 for (
size_t i = 0;
i < measurementsOnTrack->
size();
i++) {
1232 measurementsOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[pixcl->
prepRawData()->
identify()]);
1235 measurementsOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[sctcl->
prepRawData()->
identify()]);
1237 measurementsOnTrack_pixcl_sctcl_index.push_back(-1);
1241 if (outliersOnTrack) {
1242 for (
size_t i = 0;
i < outliersOnTrack->
size();
i++) {
1247 outliersOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[pixcl->
prepRawData()->
identify()]);
1249 outliersOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[sctcl->
prepRawData()->
identify()]);
1251 outliersOnTrack_pixcl_sctcl_index.push_back(-1);
1255 if (
found != trackTruthCollection->end()) {
1256 TTCindex =
found->first.index();
1257 TTCevent_index =
found->second.particleLink().eventIndex();
1258 TTCparticle_link =
found->second.particleLink().barcode();
1259 TTCprobability =
found->second.probability();
1261 TTCindex = TTCevent_index = TTCparticle_link = -999;
1262 TTCprobability = -1;
1270 (*m_TRKproperties).push_back(v_properties);
1271 (*m_TRKpattern).push_back(v_pattern);
1274 (*m_TRKmeasurementsOnTrack_pixcl_sctcl_index).push_back(measurementsOnTrack_pixcl_sctcl_index);
1275 (*m_TRKoutliersOnTrack_pixcl_sctcl_index).push_back(outliersOnTrack_pixcl_sctcl_index);
1277 (*m_TRKperigee_position).push_back(position);
1278 (*m_TRKperigee_momentum).push_back(
momentum);
1298 if (not detailedTrackTruthCollectionHandle.isValid()) {
1300 return StatusCode::FAILURE;
1302 detailedTrackTruthCollection = detailedTrackTruthCollectionHandle.cptr();
1306 (*m_DTTtrajectory_eventindex).clear();
1307 (*m_DTTtrajectory_barcode).clear();
1308 (*m_DTTstTruth_subDetType).clear();
1309 (*m_DTTstTrack_subDetType).clear();
1310 (*m_DTTstCommon_subDetType).clear();
1314 DetailedTrackTruthCollection::const_iterator detailedTrackTruthIterator = (*detailedTrackTruthCollection).begin();
1315 for (; detailedTrackTruthIterator != (*detailedTrackTruthCollection).end(); ++detailedTrackTruthIterator) {
1316 std::vector<int> DTTtrajectory_eventindex, DTTtrajectory_barcode, DTTstTruth_subDetType, DTTstTrack_subDetType,
1317 DTTstCommon_subDetType;
1318 const TruthTrajectory &traj = detailedTrackTruthIterator->second.trajectory();
1319 for (
size_t j = 0; j < traj.size(); j++) {
1320 DTTtrajectory_eventindex.push_back(traj[j].
eventIndex());
1321 DTTtrajectory_barcode.push_back(traj[j].
barcode());
1339 (*m_DTTtrajectory_eventindex).push_back(DTTtrajectory_eventindex);
1340 (*m_DTTtrajectory_barcode).push_back(DTTtrajectory_barcode);
1341 (*m_DTTstTruth_subDetType).push_back(DTTstTruth_subDetType);
1342 (*m_DTTstTrack_subDetType).push_back(DTTstTrack_subDetType);
1343 (*m_DTTstCommon_subDetType).push_back(DTTstCommon_subDetType);
1354 return StatusCode::SUCCESS;
1.8.18