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);
942 bool isRemoved(
false);
947 while ( !isRemoved && nextCellIter != myNextCells.end() ) {
948 if ( (*nextCellIter) == neighborCell ) {
949 nextCellIter = myNextCells.erase(nextCellIter);
957 pNCell->setSecondCaloTopoTmpHashCluster(myCluster);
958 nextSharedCellList.push_back(neighborCell);
959 otherCluster->remove(neighborCell);
963 if ( myCluster != otherCluster ) {
964 HashCluster *toKill = otherCluster;
965 HashCluster *toKeep = myCluster;
967 if ( !toKill->getContainsLocalMax() && toKill->size() == 1) {
977 for(;clusCellIter!=clusCellIterEnd;++clusCellIter) {
978 clusCellIter->setCaloTopoTmpHashCluster(toKeep);
980 toKeep->add(*toKill);
986 toKeep->add(neighborCell);
987 pNCell->setCaloTopoTmpHashCluster(toKeep);
993 mySeedCells.swap (myNextCells);
996 std::sort(mySeedCells.begin(),mySeedCells.end(),compareE);
999 std::sort(nextSharedCellList.begin(),nextSharedCellList.end(),compareE);
1000 if (sharedCellList.empty())
1001 sharedCellList.swap (nextSharedCellList);
1003 sharedCellList.insert(sharedCellList.end(),
1004 nextSharedCellList.begin(),
1005 nextSharedCellList.end());
1007 nextSharedCellList.clear();
1017 mySeedCells = sharedCellList;
1018 while ( !mySeedCells.empty() ) {
1021 myNextCells.clear();
1022 hashCellIter= mySeedCells.begin();
1023 hashCellIterEnd=mySeedCells.end();
1026 for(;hashCellIter!=hashCellIterEnd;++hashCellIter) {
1036 theNeighbors.clear();
1053 << (
unsigned int)hashid <<
"|"
1055 <<
"] has " << theNeighbors.size() <<
" neighbors:"
1059 for (
unsigned int iN=0;iN<theNeighbors.size();iN++) {
1066 << (
unsigned int) nId <<
"|"
1069 theNNeighbors.clear();
1071 bool foundId (
false);
1073 for (
unsigned int iNN=0;iNN<theNNeighbors.size();iNN++) {
1075 if (nOtherSubDet == ((
int)(mySubDet)) &&
1076 theNNeighbors[iNN] == hashid)
1084 msg(MSG::ERROR) <<
" Shared Cell [" << mySubDet <<
"|"
1085 << (
unsigned int)hashid <<
"|"
1087 <<
"] has bad neighbor cell[";
1090 msg() << otherSubDet <<
"|" << nId <<
"|"
1095 HashCell neighborCell = cellVector[(
unsigned int)nId -
m_hashMin];
1097 if ( pNCell && pNCell->getParentClusterIndex() == iParent && !pNCell->getShared() && !pNCell->getUsed() ) {
1102 if ( myCluster != otherCluster && mySecondCluster != otherCluster) {
1104 pNCell->setShared();
1105 myNextCells.push_back(neighborCell);
1106 sharedCellList.push_back(neighborCell);
1108 otherCluster->removeAll();
1109 pNCell->setCaloTopoTmpHashCluster(myCluster);
1110 pNCell->setSecondCaloTopoTmpHashCluster(mySecondCluster);
1115 mySeedCells.swap (myNextCells);
1118 std::sort(mySeedCells.begin(),mySeedCells.end(),compareE);
1135 hashCellIter= sharedCellList.begin();
1136 hashCellIterEnd=sharedCellList.end();
1137 for(;hashCellIter!=hashCellIterEnd;++hashCellIter) {
1146 Vector3D<double> thisPos(itrCell->
x(),itrCell->
y(),itrCell->
z());
1149 double d1 = (thisPos-
c1).
mag();
1150 double d2 = (thisPos-
c2).
mag();
1160 const double real_r_exp =
r > 10. ? 10. :
r < -10. ? -10. :
r;
1161 const double real_r =
exp(real_r_exp);
1162 const double r_reverse =
exp(-real_r_exp);
1164 const double reverse_weight =
e2 / (
e2 +
e1 * r_reverse);
1188 hashCellIter= sharedCellList.begin();
1189 hashCellIterEnd=sharedCellList.end();
1190 for(;hashCellIter!=hashCellIterEnd;++hashCellIter) {
1194 firstCluster->
add(*hashCellIter);
1195 secondCluster->add(*hashCellIter);
1197 nShared = sharedCellList.size();
1204 std::vector<std::unique_ptr<CaloProtoCluster> > myCaloClusters;
1205 myCaloClusters.reserve (500);
1206 std::vector<std::unique_ptr<CaloProtoCluster> > myRestClusters;
1207 myRestClusters.resize(clusColl->size());
1210 for (;hashClusIter!=hashClusIterEnd;++hashClusIter) {
1211 HashCluster * tmpCluster = (*hashClusIter);
1212 if ( (
m_useGPUCriteria && tmpCluster->getContainsLocalMax()) || tmpCluster->size() > 1 ) {
1227 std::unique_ptr<CaloProtoCluster> myCluster = std::make_unique<CaloProtoCluster>(myCellCollLink);
1228 ATH_MSG_DEBUG(
"[CaloCluster@" << myCluster.get() <<
"] created in <myCaloClusters>.");
1232 for(;clusCellIter!=clusCellIterEnd;++clusCellIter) {
1235 double myWeight = itrCell.
weight();
1239 if (shared_weight < 0.) {
1241 myWeight *= 1.0 + shared_weight;
1243 myWeight *= -shared_weight;
1247 myWeight *= shared_weight;
1249 myWeight *= 1.0 - shared_weight;
1261 ATH_MSG_DEBUG(
"[CaloCluster@" << myCluster.get() <<
"] size: " << myCluster->
size());
1262 myCaloClusters.push_back(std::move(myCluster));
1264 else if ( tmpCluster->size() == 1 ) {
1269 if ( hasLocalMaxVector[tmpCluster->getParentClusterIndex()]) {
1273 if (!myRestClusters[tmpCluster->getParentClusterIndex()]) {
1274 myRestClusters[tmpCluster->getParentClusterIndex()] = std::make_unique<CaloProtoCluster>(myCellColl);
1276 ATH_MSG_DEBUG(
"[CaloCluster@" << myRestClusters[tmpCluster->getParentClusterIndex()].get()
1277 <<
"] created in <myRestClusters>");
1278 myRestClusters[tmpCluster->getParentClusterIndex()]->getCellLinks()->reserve(tmpCluster->size());
1281 for(;clusCellIter!=clusCellIterEnd;++clusCellIter) {
1284 const double myWeight = itrCell.
weight();
1285 myRestClusters[tmpCluster->getParentClusterIndex()]->addCell(itrCell.
index(),myWeight);
1287 ATH_MSG_DEBUG(
"[CaloCluster@" << myRestClusters[tmpCluster->getParentClusterIndex()].get()
1288 <<
"] size: " << myRestClusters[tmpCluster->getParentClusterIndex()]->size());
1296 clusCollIter = clusColl->begin();
1297 for (; clusCollIter != clusCollIterEnd; ++clusCollIter,++iClusterNumber){
1299 if ( !hasLocalMaxVector[iClusterNumber] ) {
1301 myCaloClusters.push_back(std::make_unique<CaloProtoCluster>(parentCluster->
getCellLinks()));
1302 ATH_MSG_DEBUG(
"[CaloProtoCluster@" << myCaloClusters.back().get() <<
"] with "
1303 << myCaloClusters.back()->size() <<
"cells cloned from "
1304 << parentCluster <<
" with " << parentCluster->
size() <<
" cells");
1306 else if (myRestClusters[iClusterNumber]) {
1308 <<
"] pushed into <myCaloClusters> with "
1309 << myRestClusters[iClusterNumber]->
size() <<
" cells");
1310 myCaloClusters.push_back(std::move(myRestClusters[iClusterNumber]));
1314 std::sort(myCaloClusters.begin(),myCaloClusters.end(),[](
const std::unique_ptr<CaloProtoCluster>& pc1,
1315 const std::unique_ptr<CaloProtoCluster>& pc2) {
1318 volatile double et1(pc1->et());
1319 volatile double et2(pc2->et());
1327 clusColl->reserve (myCaloClusters.size());
1333 for(
const auto& protoCluster : myCaloClusters) {
1335 clusColl->push_back(xAODCluster);
1336 xAODCluster->
addCellLink(protoCluster->releaseCellLinks());
1339 ATH_MSG_DEBUG(
"CaloCluster@" << xAODCluster <<
" pushed into "
1340 <<
"CaloClusterContainer@" << clusColl);
1343 <<
"->size() = " << clusColl->size());
1345 <<
" MeV, Et = " << xAODCluster->
et()
1346 <<
" MeV, NCells = " << xAODCluster->
size());
1347 eTot+=xAODCluster->
e();
1348 nTot+=xAODCluster->
size();
1350 if ( std::abs(xAODCluster->
e()) > eMax )
1351 eMax = std::abs(xAODCluster->
e());
1354 <<
" MeV, NCells = " << nTot
1355 <<
" (including NShared = " << nShared <<
" twice)");
1357 if ( std::abs(
eTot) > eMax )
1358 eMax = std::abs(
eTot);
1359 if ( std::abs(
eTot-eTotOrig)>0.001*eMax ){
1360 msg(MSG::WARNING) <<
"Energy sum for split Clusters = " <<
eTot <<
" MeV does not equal original sum = " << eTotOrig <<
" MeV !" <<
endmsg;
1362 if ( abs(nTot-nShared-nTotOrig) > 0 ) {
1363 msg(MSG::ERROR) <<
"Cell sum for split Clusters does not equal original sum!" <<
endmsg;
1366 tmpcell_pool.
erase();
1367 tmpclus_pool.
erase();
1369 return StatusCode::SUCCESS;
1.8.18