Execute on an entire collection of clusters.
324 std::vector<HashCell> mySeedCells;
325 mySeedCells.reserve (200);
327 std::vector<HashCluster *> myHashClusters;
328 myHashClusters.reserve (20000);
330 std::vector<CaloTopoSplitterClusterCell *> allCellList;
331 allCellList.reserve (20000);
338 std::vector<int> hasLocalMaxVector;
339 hasLocalMaxVector.resize(clusColl->size(),0);
340 int iClusterNumber = 0;
349 if (clusCollIter != clusCollIterEnd) {
350 clusterSize = (*clusCollIter)->clusterSize();
354 msg(MSG::ERROR) <<
"Can't get valid links to CaloCells (CaloClusterCellLink)!" <<
endmsg;
355 return StatusCode::FAILURE;
361 return StatusCode::SUCCESS;
365 for (; clusCollIter != clusCollIterEnd; ++clusCollIter, ++iClusterNumber ){
369 msg(MSG::ERROR) <<
"Can't get valid links to CaloCells (CaloClusterCellLink)!" <<
endmsg;
370 return StatusCode::FAILURE;
373 eTotOrig+=parentCluster->
e();
374 nTotOrig+=cellLinks->
size();
384 for(; cellIter != cellIterEnd; cellIter++ ){
387 const CaloCell* pCell = (*cellIter);
392 float signedRatio = 0;
393 bool is_secondary =
false;
405 signedRatio = pCell->
e();
415 signedRatio = -pCell->
e();
425 size_t myIndex=cellIter.
index();
430 signedRatio,parentCluster,
431 iClusterNumber,is_secondary);
435 msg(MSG::INFO) <<
" [ExtId|Id|SubDet|HashId|eta|phi|iParent|E]: "
439 <<
"|" << (
unsigned int)hashid
440 <<
"|" << pCell->
eta()
441 <<
"|" << pCell->
phi()
442 <<
"|" << iClusterNumber
443 <<
"|" << signedRatio
446 HashCell hashCell(tmpClusterCell);
447 HashCluster *tmpCluster =
448 new (tmpclus_pool.
allocate()) HashCluster (tmplist_pool);
450 tmpCluster->add(hashCell);
451 myHashClusters.push_back(tmpCluster);
452 allCellList.push_back(tmpClusterCell);
453 cellVector[(
unsigned int)hashid -
m_hashMin] = hashCell;
460 std::vector<IdentifierHash> theNeighbors;
461 theNeighbors.reserve(22);
462 std::vector<IdentifierHash> theSuperNeighbors;
463 theSuperNeighbors.reserve(22);
464 std::vector<IdentifierHash> theNNeighbors;
465 theNNeighbors.reserve(22);
466 std::vector<IdentifierHash> theCurrentNeighbors;
467 theCurrentNeighbors.reserve(88);
468 std::vector<IdentifierHash> theNextNeighbors;
469 theNextNeighbors.reserve(88);
474 for(;allCellIter != allCellIterEnd;++allCellIter) {
480 if(
m_absOpt) myEnergy=std::abs(myEnergy);
486 theNeighbors.clear();
488 for (
unsigned int iN=0;iN<theNeighbors.size();iN++) {
490 HashCell neighborCell = cellVector[(
unsigned int)nId -
m_hashMin];
492 if ( pNeighCell && pNeighCell->getParentClusterIndex() == iParent) {
494 if ( ((myEnergy > pNeighCell->getSignedRatio() ) && !
m_absOpt) ||
495 (
m_absOpt && myEnergy > std::abs(pNeighCell->getSignedRatio() ) ) ||
496 pNeighCell->getSecondary() ) {
500 if (std::abs(pNeighCell->getSignedRatio()) < myEnergy )
501 pNeighCell->setUsed();
511 mySeedCells.push_back(cellVector[(
unsigned int)hashid -
m_hashMin]);
512 hasLocalMaxVector[iParent]++;
518 std::vector<CaloTopoSplitterClusterCell *> myPotentialSecondarySeeds;
524 myPotentialSecondarySeeds.reserve(100);
528 allCellIter=allCellList.begin();
529 for(;allCellIter != allCellIterEnd;++allCellIter) {
535 if(
m_absOpt) myEnergy=std::abs(myEnergy);
542 theNeighbors.clear();
544 for (
unsigned int iN=0;iN<theNeighbors.size();iN++) {
546 HashCell neighborCell = cellVector[(
unsigned int)nId -
m_hashMin];
548 if ( pNeighCell && pNeighCell->getParentClusterIndex() == iParent) {
550 if ( std::abs(myEnergy) > std::abs(pNeighCell->getSignedRatio()) ) {
551 pNeighCell->setUsed();
563 myPotentialSecondarySeeds.push_back(pClusCell);
576 theCurrentNeighbors.clear();
577 theCurrentNeighbors.push_back(hashid);
578 while (
isLocalMax && !theCurrentNeighbors.empty() ) {
581 theNextNeighbors.clear();
582 for (
unsigned int iN=0;iN<theCurrentNeighbors.size()
584 theNeighbors.clear();
585 theSuperNeighbors.clear();
590 theNeighbors.insert(theNeighbors.end(),theSuperNeighbors.begin(),theSuperNeighbors.end());
591 for(
unsigned int iNN=0;iNN<theNeighbors.size() &&
isLocalMax;iNN++) {
596 hashCellIter= mySeedCells.begin();
597 hashCellIterEnd=mySeedCells.end();
600 for(;hashCellIter!=hashCellIterEnd
602 if ( cellVector[(
unsigned int)nId -
m_hashMin]
612 bool doInclude(
true);
613 for(
unsigned int iNNN=0;iNNN<theNextNeighbors.size();iNNN++) {
614 if ( theNextNeighbors[iNNN] == theNeighbors[iNN] ) {
620 theNextNeighbors.push_back(theNeighbors[iNN]);
624 theCurrentNeighbors.swap (theNextNeighbors);
631 std::vector<IdentifierHash> theCurrentNeighbors;
632 std::vector<IdentifierHash> theNextNeighbors;
634 theCurrentNeighbors.push_back(hashid);
635 while (
isLocalMax && !theCurrentNeighbors.empty() ) {
638 theNextNeighbors.clear();
639 for (
unsigned int iN=0;iN<theCurrentNeighbors.size()
641 theNeighbors.clear();
642 theSuperNeighbors.clear();
647 theNeighbors.insert(theNeighbors.end(),theSuperNeighbors.begin(),theSuperNeighbors.end());
648 for(
unsigned int iNN=0;iNN<theNeighbors.size() &&
isLocalMax;iNN++) {
653 hashCellIter= mySeedCells.begin();
654 hashCellIterEnd=mySeedCells.end();
657 for(;hashCellIter!=hashCellIterEnd
659 if ( cellVector[(
unsigned int)nId -
m_hashMin]
669 bool doInclude(
true);
670 for(
unsigned int iNNN=0;iNNN<theNextNeighbors.size();iNNN++) {
671 if ( theNextNeighbors[iNNN] == theNeighbors[iNN] ) {
677 theNextNeighbors.push_back(theNeighbors[iNN]);
681 theCurrentNeighbors.swap (theNextNeighbors);
687 mySeedCells.push_back(cellVector[(
unsigned int)hashid -
m_hashMin]);
688 hasLocalMaxVector[iParent]++;
694 allCellIter=allCellList.begin();
695 for(;allCellIter != allCellIterEnd;++allCellIter) {
696 (*allCellIter)->setUnused();
698 (*allCellIter)->setSignedRatio(itrCell->
e());
707 std::sort(myPotentialSecondarySeeds.begin(), myPotentialSecondarySeeds.end(), [&](
const auto &
a,
const auto &
b) {
708 return compareEWithIndex(cellVector[(unsigned int)a->getID() - m_hashMin], cellVector[(unsigned int)b->getID() - m_hashMin]);
715 std::vector<bool> secondarySeedExclude(myPotentialSecondarySeeds.size(),
false);
717 for (
unsigned int i = 0;
i < myPotentialSecondarySeeds.size(); ++
i) {
727 theCurrentNeighbors.clear();
728 theCurrentNeighbors.push_back(hashid);
730 while (
isLocalMax && !theCurrentNeighbors.empty() ) {
732 theNextNeighbors.clear();
733 for (
const IdentifierHash & currentNeighbor : theCurrentNeighbors) {
735 theNeighbors.clear();
736 theSuperNeighbors.clear();
746 theNeighbors.insert(theNeighbors.end(), theSuperNeighbors.begin(), theSuperNeighbors.end());
750 for (
const auto & seedCell: mySeedCells) {
751 if (cellVector[(
unsigned int)nId -
m_hashMin] == seedCell) {
756 for (
unsigned int j = 0; j <
i; ++j) {
757 if (cellVector[(
unsigned int)nId -
m_hashMin] == myPotentialSecondarySeeds[j]) {
762 bool doInclude(
true);
765 if (nextNId == nId) {
772 theNextNeighbors.push_back(nId);
777 theCurrentNeighbors.swap (theNextNeighbors);
786 secondarySeedExclude[
i] =
true;
790 for (
unsigned int i = 0;
i < myPotentialSecondarySeeds.size(); ++
i) {
794 if ( !secondarySeedExclude[
i] ) {
795 mySeedCells.push_back(cellVector[(
unsigned int)hashid -
m_hashMin]);
796 hasLocalMaxVector[iParent]++;
802 std::vector<HashCell> sharedCellList;
803 std::vector<HashCell> nextSharedCellList;
808 hashCellIter= mySeedCells.begin();
809 hashCellIterEnd=mySeedCells.end();
812 for(;hashCellIter!=hashCellIterEnd;++hashCellIter) {
813 hashCellIter->getCaloTopoTmpClusterCell()->setUsed();
814 HashCluster *myCluster = hashCellIter->getCaloTopoTmpClusterCell()->getCaloTopoTmpHashCluster();
815 myCluster->setContainsLocalMax();
824 auto compareE = [&,
this](
auto && ... ps) {
826 return compareEWithIndex(std::forward<decltype(ps)>(ps)...);
829 return compareEOriginal(std::forward<decltype(ps)>(ps)...);
833 std::sort(mySeedCells.begin(),mySeedCells.end(),compareE);
836 hashCellIter= mySeedCells.begin();
837 hashCellIterEnd=mySeedCells.end();
840 for(;hashCellIter!=hashCellIterEnd;++hashCellIter) {
842 << hashCellIter->getCaloTopoTmpClusterCell()->getSubDet()
844 << (
unsigned int)hashCellIter->getCaloTopoTmpClusterCell()->getID()
846 << hashCellIter->getCaloTopoTmpClusterCell()->getSignedRatio()
851 std::vector<HashCell> myNextCells;
852 myNextCells.reserve (4096);
853 while ( !mySeedCells.empty() ) {
856 hashCellIter= mySeedCells.begin();
857 hashCellIterEnd=mySeedCells.end();
860 for(;hashCellIter!=hashCellIterEnd;++hashCellIter) {
869 theNeighbors.clear();
886 << (
unsigned int)hashid <<
"|"
888 <<
"] has " << theNeighbors.size() <<
" neighbors:"
892 for (
unsigned int iN=0;iN<theNeighbors.size();iN++) {
899 << (
unsigned int) nId <<
"|"
901 theNNeighbors.clear();
903 bool foundId (
false);
905 for (
unsigned int iNN=0;iNN<theNNeighbors.size();iNN++) {
907 if (nOtherSubDet == ((
int)(mySubDet)) &&
908 theNNeighbors[iNN] == hashid) {
915 msg(MSG::ERROR) <<
" Cell [" << mySubDet <<
"|"
916 << (
unsigned int)hashid <<
"|"
918 <<
"] has bad neighbor cell[";
921 msg() << otherSubDet <<
"|" << nId <<
"|"
926 HashCell neighborCell = cellVector[(
unsigned int)nId -
m_hashMin];
928 if ( pNCell && pNCell->getParentClusterIndex() == iParent && !pNCell->getShared() ) {
931 if ( !pNCell->getUsed() ) {
933 myNextCells.push_back(neighborCell);
943 bool isRemoved(
false);
949 while ( !isRemoved && nextCellIter != nextCellIterEnd ) {
950 if ( (*nextCellIter) == neighborCell ) {
951 nextCellIter = myNextCells.erase(nextCellIter);
959 pNCell->setSecondCaloTopoTmpHashCluster(myCluster);
960 nextSharedCellList.push_back(neighborCell);
961 otherCluster->remove(neighborCell);
965 if ( myCluster != otherCluster ) {
966 HashCluster *toKill = otherCluster;
967 HashCluster *toKeep = myCluster;
969 if ( !toKill->getContainsLocalMax() && toKill->size() == 1) {
979 for(;clusCellIter!=clusCellIterEnd;++clusCellIter) {
980 clusCellIter->setCaloTopoTmpHashCluster(toKeep);
982 toKeep->add(*toKill);
988 toKeep->add(neighborCell);
989 pNCell->setCaloTopoTmpHashCluster(toKeep);
995 mySeedCells.swap (myNextCells);
998 std::sort(mySeedCells.begin(),mySeedCells.end(),compareE);
1001 std::sort(nextSharedCellList.begin(),nextSharedCellList.end(),compareE);
1002 if (sharedCellList.empty())
1003 sharedCellList.swap (nextSharedCellList);
1005 sharedCellList.insert(sharedCellList.end(),
1006 nextSharedCellList.begin(),
1007 nextSharedCellList.end());
1009 nextSharedCellList.clear();
1019 mySeedCells = sharedCellList;
1020 while ( !mySeedCells.empty() ) {
1023 myNextCells.clear();
1024 hashCellIter= mySeedCells.begin();
1025 hashCellIterEnd=mySeedCells.end();
1028 for(;hashCellIter!=hashCellIterEnd;++hashCellIter) {
1038 theNeighbors.clear();
1055 << (
unsigned int)hashid <<
"|"
1057 <<
"] has " << theNeighbors.size() <<
" neighbors:"
1061 for (
unsigned int iN=0;iN<theNeighbors.size();iN++) {
1068 << (
unsigned int) nId <<
"|"
1071 theNNeighbors.clear();
1073 bool foundId (
false);
1075 for (
unsigned int iNN=0;iNN<theNNeighbors.size();iNN++) {
1077 if (nOtherSubDet == ((
int)(mySubDet)) &&
1078 theNNeighbors[iNN] == hashid)
1086 msg(MSG::ERROR) <<
" Shared Cell [" << mySubDet <<
"|"
1087 << (
unsigned int)hashid <<
"|"
1089 <<
"] has bad neighbor cell[";
1092 msg() << otherSubDet <<
"|" << nId <<
"|"
1097 HashCell neighborCell = cellVector[(
unsigned int)nId -
m_hashMin];
1099 if ( pNCell && pNCell->getParentClusterIndex() == iParent && !pNCell->getShared() && !pNCell->getUsed() ) {
1104 if ( myCluster != otherCluster && mySecondCluster != otherCluster) {
1106 pNCell->setShared();
1107 myNextCells.push_back(neighborCell);
1108 sharedCellList.push_back(neighborCell);
1110 otherCluster->removeAll();
1111 pNCell->setCaloTopoTmpHashCluster(myCluster);
1112 pNCell->setSecondCaloTopoTmpHashCluster(mySecondCluster);
1117 mySeedCells.swap (myNextCells);
1120 std::sort(mySeedCells.begin(),mySeedCells.end(),compareE);
1137 hashCellIter= sharedCellList.begin();
1138 hashCellIterEnd=sharedCellList.end();
1139 for(;hashCellIter!=hashCellIterEnd;++hashCellIter) {
1148 Vector3D<double> thisPos(itrCell->
x(),itrCell->
y(),itrCell->
z());
1151 double d1 = (thisPos-
c1).
mag();
1152 double d2 = (thisPos-
c2).
mag();
1162 const double real_r_exp =
r > 10. ? 10. :
r < -10. ? -10. :
r;
1163 const double real_r =
exp(real_r_exp);
1164 const double r_reverse =
exp(-real_r_exp);
1166 const double reverse_weight =
e2 / (
e2 +
e1 * r_reverse);
1190 hashCellIter= sharedCellList.begin();
1191 hashCellIterEnd=sharedCellList.end();
1192 for(;hashCellIter!=hashCellIterEnd;++hashCellIter) {
1196 firstCluster->
add(*hashCellIter);
1197 secondCluster->add(*hashCellIter);
1199 nShared = sharedCellList.size();
1206 std::vector<std::unique_ptr<CaloProtoCluster> > myCaloClusters;
1207 myCaloClusters.reserve (500);
1208 std::vector<std::unique_ptr<CaloProtoCluster> > myRestClusters;
1209 myRestClusters.resize(clusColl->size());
1212 for (;hashClusIter!=hashClusIterEnd;++hashClusIter) {
1213 HashCluster * tmpCluster = (*hashClusIter);
1214 if ( (
m_useGPUCriteria && tmpCluster->getContainsLocalMax()) || tmpCluster->size() > 1 ) {
1229 std::unique_ptr<CaloProtoCluster> myCluster = std::make_unique<CaloProtoCluster>(myCellCollLink);
1230 ATH_MSG_DEBUG(
"[CaloCluster@" << myCluster.get() <<
"] created in <myCaloClusters>.");
1234 for(;clusCellIter!=clusCellIterEnd;++clusCellIter) {
1237 double myWeight = itrCell.
weight();
1241 if (shared_weight < 0.) {
1243 myWeight *= 1.0 + shared_weight;
1245 myWeight *= -shared_weight;
1249 myWeight *= shared_weight;
1251 myWeight *= 1.0 - shared_weight;
1263 ATH_MSG_DEBUG(
"[CaloCluster@" << myCluster.get() <<
"] size: " << myCluster->
size());
1264 myCaloClusters.push_back(std::move(myCluster));
1266 else if ( tmpCluster->size() == 1 ) {
1271 if ( hasLocalMaxVector[tmpCluster->getParentClusterIndex()]) {
1275 if (!myRestClusters[tmpCluster->getParentClusterIndex()]) {
1276 myRestClusters[tmpCluster->getParentClusterIndex()] = std::make_unique<CaloProtoCluster>(myCellColl);
1278 ATH_MSG_DEBUG(
"[CaloCluster@" << myRestClusters[tmpCluster->getParentClusterIndex()].get()
1279 <<
"] created in <myRestClusters>");
1280 myRestClusters[tmpCluster->getParentClusterIndex()]->getCellLinks()->reserve(tmpCluster->size());
1283 for(;clusCellIter!=clusCellIterEnd;++clusCellIter) {
1286 const double myWeight = itrCell.
weight();
1287 myRestClusters[tmpCluster->getParentClusterIndex()]->addCell(itrCell.
index(),myWeight);
1289 ATH_MSG_DEBUG(
"[CaloCluster@" << myRestClusters[tmpCluster->getParentClusterIndex()].get()
1290 <<
"] size: " << myRestClusters[tmpCluster->getParentClusterIndex()]->size());
1298 clusCollIter = clusColl->begin();
1299 for (; clusCollIter != clusCollIterEnd; ++clusCollIter,++iClusterNumber){
1301 if ( !hasLocalMaxVector[iClusterNumber] ) {
1303 myCaloClusters.push_back(std::make_unique<CaloProtoCluster>(parentCluster->
getCellLinks()));
1304 ATH_MSG_DEBUG(
"[CaloProtoCluster@" << myCaloClusters.back().get() <<
"] with "
1305 << myCaloClusters.back()->size() <<
"cells cloned from "
1306 << parentCluster <<
" with " << parentCluster->
size() <<
" cells");
1308 else if (myRestClusters[iClusterNumber]) {
1310 <<
"] pushed into <myCaloClusters> with "
1311 << myRestClusters[iClusterNumber]->
size() <<
" cells");
1312 myCaloClusters.push_back(std::move(myRestClusters[iClusterNumber]));
1316 std::sort(myCaloClusters.begin(),myCaloClusters.end(),[](
const std::unique_ptr<CaloProtoCluster>& pc1,
1317 const std::unique_ptr<CaloProtoCluster>& pc2) {
1320 volatile double et1(pc1->et());
1321 volatile double et2(pc2->et());
1329 clusColl->reserve (myCaloClusters.size());
1335 for(
const auto& protoCluster : myCaloClusters) {
1337 clusColl->push_back(xAODCluster);
1338 xAODCluster->
addCellLink(protoCluster->releaseCellLinks());
1341 ATH_MSG_DEBUG(
"CaloCluster@" << xAODCluster <<
" pushed into "
1342 <<
"CaloClusterContainer@" << clusColl);
1345 <<
"->size() = " << clusColl->size());
1347 <<
" MeV, Et = " << xAODCluster->
et()
1348 <<
" MeV, NCells = " << xAODCluster->
size());
1349 eTot+=xAODCluster->
e();
1350 nTot+=xAODCluster->
size();
1352 if ( std::abs(xAODCluster->
e()) > eMax )
1353 eMax = std::abs(xAODCluster->
e());
1356 <<
" MeV, NCells = " << nTot
1357 <<
" (including NShared = " << nShared <<
" twice)");
1359 if ( std::abs(
eTot) > eMax )
1360 eMax = std::abs(
eTot);
1361 if ( std::abs(
eTot-eTotOrig)>0.001*eMax ){
1362 msg(MSG::WARNING) <<
"Energy sum for split Clusters = " <<
eTot <<
" MeV does not equal original sum = " << eTotOrig <<
" MeV !" <<
endmsg;
1364 if ( abs(nTot-nShared-nTotOrig) > 0 ) {
1365 msg(MSG::ERROR) <<
"Cell sum for split Clusters does not equal original sum!" <<
endmsg;
1368 tmpcell_pool.
erase();
1369 tmpclus_pool.
erase();
1371 return StatusCode::SUCCESS;
1.8.18