diff --git a/RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h b/RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h index d9cf9ffce395d..3b8b5db320b08 100644 --- a/RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h +++ b/RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h @@ -20,6 +20,10 @@ #include "FWCore/Framework/interface/ConsumesCollector.h" #include "PhysicsTools/TensorFlow/interface/TensorFlow.h" #include "CommonTools/RecoAlgos/interface/MultiVectorManager.h" +#include "MagneticField/Engine/interface/MagneticField.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" +#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h" namespace edm { class Event; @@ -53,6 +57,11 @@ namespace ticl { std::vector>& linkedResultTracksters, std::vector>& linkedTracksterIdToInputTracksterId) = 0; + virtual void initialize(const HGCalDDDConstants* hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) = 0; + static void fillPSetDescription(edm::ParameterSetDescription& desc) { desc.add("algo_verbosity", 0); }; protected: diff --git a/RecoHGCal/TICL/plugins/TICLGraph.cc b/RecoHGCal/TICL/plugins/TICLGraph.cc new file mode 100644 index 0000000000000..6df936589804a --- /dev/null +++ b/RecoHGCal/TICL/plugins/TICLGraph.cc @@ -0,0 +1,57 @@ +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "TICLGraph.h" + +void Node::findSubComponents(std::vector& graph, std::vector& subComponent, std::string tabs) { + tabs += "\t"; + if (!alreadyVisited_) { + LogDebug("TICLGraph") << tabs << " Visiting node " << index_ << std::endl; + alreadyVisited_ = true; + subComponent.push_back(index_); + for (auto const& neighbour : neighboursId_) { + LogDebug("TICLGraph") << tabs << " Trying to visit " << neighbour << std::endl; + graph[neighbour].findSubComponents(graph, subComponent, tabs); + } + } +} + +std::vector> TICLGraph::findSubComponents() { + std::vector> components; + for (auto const& node : nodes_) { + auto const id = node.getId(); + if (isRootNode_[id]) { + //LogDebug("TICLGraph") << "DFS Starting From " << id << std::endl; + std::string tabs = "\t"; + std::vector tmpSubComponents; + nodes_[id].findSubComponents(nodes_, tmpSubComponents, tabs); + components.push_back(tmpSubComponents); + } + } + return components; +} + +void TICLGraph::dfsForCC(unsigned int nodeIndex, + std::unordered_set& visited, + std::vector& component) const { + visited.insert(nodeIndex); + component.push_back(nodeIndex); + + for (auto const& neighbourIndex : nodes_[nodeIndex].getNeighbours()) { + if (visited.find(neighbourIndex) == visited.end()) { + dfsForCC(neighbourIndex, visited, component); + } + } +} + +std::vector> TICLGraph::getConnectedComponents() const { + std::unordered_set visited; + std::vector> components; + + for (unsigned int i = 0; i < nodes_.size(); ++i) { + if (visited.find(i) == visited.end()) { + std::vector component; + dfsForCC(i, visited, component); + components.push_back(component); + } + } + return components; +} diff --git a/RecoHGCal/TICL/plugins/TICLGraph.h b/RecoHGCal/TICL/plugins/TICLGraph.h new file mode 100644 index 0000000000000..5956c0044b29c --- /dev/null +++ b/RecoHGCal/TICL/plugins/TICLGraph.h @@ -0,0 +1,58 @@ +#ifndef DataFormats_HGCalReco_TICLGraph_h +#define DataFormats_HGCalReco_TICLGraph_h + +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include + +class Node { +public: + Node() = default; + Node(unsigned index, bool isTrackster = true) : index_(index), isTrackster_(isTrackster), alreadyVisited_{false} {}; + + inline void addNeighbour(unsigned int trackster_id) { neighboursId_.push_back(trackster_id); } + + inline const unsigned int getId() const { return index_; } + std::vector getNeighbours() const { return neighboursId_; } + void findSubComponents(std::vector& graph, std::vector& subComponent, std::string tabs); + + ~Node() = default; + +private: + unsigned index_; + bool isTrackster_; + + std::vector neighboursId_; + bool alreadyVisited_; + + //bool areCompatible(const std::vector& graph, const unsigned int& outerNode) { return true; }; +}; + +class TICLGraph { +public: + // can i remove default constructor ?? edm::Wrapper problem + // without default constructor i could initialize connectedComponents when building the Graph + TICLGraph() = default; + TICLGraph(std::vector& n, std::vector isRootNode) { + nodes_ = n; + isRootNode_ = isRootNode; + }; + inline const std::vector& getNodes() const { return nodes_; } + inline const Node& getNode(int i) const { return nodes_[i]; } + + std::vector> findSubComponents(); + + ~TICLGraph() = default; + + void dfsForCC(unsigned int nodeIndex, + std::unordered_set& visited, + std::vector& component) const; + + std::vector> getConnectedComponents() const; + +private: + std::vector nodes_; + std::vector isRootNode_; +}; + +#endif diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.cc b/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.cc index 313887b778ce1..88f7e86694bee 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.cc +++ b/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.cc @@ -1,11 +1,11 @@ -// #include "TracksterLinkingbySkeletons.h" // #include "TracksterLinkingbySuperClustering.h" #include "FWCore/ParameterSet/interface/ValidatedPluginFactoryMacros.h" #include "FWCore/ParameterSet/interface/ValidatedPluginMacros.h" #include "TracksterLinkingbyFastJet.h" +#include "TracksterLinkingbySkeletons.h" #include "RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.h" EDM_REGISTER_VALIDATED_PLUGINFACTORY(TracksterLinkingPluginFactory, "TracksterLinkingPluginFactory"); -// DEFINE_EDM_VALIDATED_PLUGIN(TracksterLinkingPluginFactory, ticl::TracksterLinkingbySkeletons, "Skeletons"); +DEFINE_EDM_VALIDATED_PLUGIN(TracksterLinkingPluginFactory, ticl::TracksterLinkingbySkeletons, "Skeletons"); // DEFINE_EDM_VALIDATED_PLUGIN(TracksterLinkingPluginFactory, ticl::TracksterLinkingbySuperClustering, "SuperClustering"); DEFINE_EDM_VALIDATED_PLUGIN(TracksterLinkingPluginFactory, ticl::TracksterLinkingbyFastJet, "FastJet"); diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.cc b/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.cc index cd92e33947bb4..0c7bf5bc0e869 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.cc +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.cc @@ -45,4 +45,4 @@ void TracksterLinkingbyFastJet::linkTracksters( linkedResultTracksters.push_back(linkedTracksters); } } -} \ No newline at end of file +} diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.h b/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.h index 3e6014d104c51..023058bfb3128 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.h +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.h @@ -37,6 +37,12 @@ namespace ticl { std::vector& resultTracksters, std::vector>& linkedResultTracksters, std::vector>& linkedTracksterIdToInputTracksterId) override; + + void initialize(const HGCalDDDConstants* hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) override{}; + static void fillPSetDescription(edm::ParameterSetDescription& iDesc) { iDesc.add("algo_verbosity", 0); iDesc.add("jet_algorithm", 2) diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.cc b/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.cc new file mode 100644 index 0000000000000..064b7bcbef051 --- /dev/null +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.cc @@ -0,0 +1,378 @@ +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/HGCalReco/interface/Common.h" +#include "DataFormats/GeometrySurface/interface/BoundDisk.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" +#include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h" +#include "RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h" +#include "RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.h" +#include "TICLGraph.h" + +using namespace ticl; + +TracksterLinkingbySkeletons::TracksterLinkingbySkeletons(const edm::ParameterSet &conf, edm::ConsumesCollector iC) + : TracksterLinkingAlgoBase(conf, iC), + timing_quality_threshold_(conf.getParameter("track_time_quality_threshold")), + del_(conf.getParameter("wind")), + angle_first_cone_(conf.getParameter("angle0")), + angle_second_cone_(conf.getParameter("angle1")), + angle_third_cone_(conf.getParameter("angle2")), + pcaQ_(conf.getParameter("pcaQuality")), + pcaQLCSize_(conf.getParameter("pcaQualityLCSize")), + dotCut_(conf.getParameter("dotProdCut")), + maxDistSkeletonsSq_(conf.getParameter("maxDistSkeletonsSq")), + max_height_cone_(conf.getParameter("maxConeHeight")) {} + +void TracksterLinkingbySkeletons::buildLayers() { + // build disks at HGCal front & EM-Had interface for track propagation + + float zVal = hgcons_->waferZ(1, true); + std::pair rMinMax = hgcons_->rangeR(zVal, true); + + float zVal_interface = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(); + std::pair rMinMax_interface = hgcons_->rangeR(zVal_interface, true); + + for (int iSide = 0; iSide < 2; ++iSide) { + float zSide = (iSide == 0) ? (-1. * zVal) : zVal; + firstDisk_[iSide] = + std::make_unique(Disk::build(Disk::PositionType(0, 0, zSide), + Disk::RotationType(), + SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5)) + .get()); + + zSide = (iSide == 0) ? (-1. * zVal_interface) : zVal_interface; + interfaceDisk_[iSide] = std::make_unique( + Disk::build(Disk::PositionType(0, 0, zSide), + Disk::RotationType(), + SimpleDiskBounds(rMinMax_interface.first, rMinMax_interface.second, zSide - 0.5, zSide + 0.5)) + .get()); + } +} + +void TracksterLinkingbySkeletons::initialize(const HGCalDDDConstants *hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) { + hgcons_ = hgcons; + rhtools_ = rhtools; + buildLayers(); + + bfield_ = bfieldH; + propagator_ = propH; +} + +float TracksterLinkingbySkeletons::findSkeletonPoints(float percentage, + const float trackster_energy, + const std::vector vertices, + const hgcal::RecHitTools &rhtools, + const std::vector &layerClusters) { + std::vector energyInLayer(rhtools_.lastLayer(), 0.); + std::vector cumulativeEnergyInLayer(rhtools_.lastLayer(), 0.); + for (auto const &v : vertices) { + auto const &lc = layerClusters[v]; + auto const &n_lay = rhtools_.getLayerWithOffset(lc.hitsAndFractions()[0].first); + energyInLayer[n_lay] += lc.energy() / trackster_energy; + } + if (TracksterLinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced) { + for (size_t iLay = 0; iLay < energyInLayer.size(); ++iLay) { + LogDebug("TracksterLinkingbySkeletons") + << "Layer " << iLay << " contains a Trackster energy fraction of " << energyInLayer[iLay] + << " and the trackster energy is " << trackster_energy << "\n"; + } + } + auto sum = 0.; + for (size_t iC = 0; iC != energyInLayer.size(); iC++) { + sum += energyInLayer[iC]; + cumulativeEnergyInLayer[iC] = sum; + } + if (TracksterLinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced) { + for (size_t iLay = 0; iLay < cumulativeEnergyInLayer.size(); ++iLay) { + LogDebug("TracksterLinkingbySkeletons") + << "Layer " << iLay << " cumulative has a Trackster energy fraction of " << cumulativeEnergyInLayer[iLay] + << " and the trackster energy is " << trackster_energy << "\n"; + } + } + auto layerI = + std::min_element(cumulativeEnergyInLayer.begin(), cumulativeEnergyInLayer.end(), [percentage](float a, float b) { + // Check if 'a' and 'b' are both greater than 0 + if (a > 0 && b > 0) { + // Compare based on absolute difference from 'percentage' + return std::abs(a - percentage) < std::abs(b - percentage); + } else if (a > 0) { + // 'a' is greater than 0, so it is the better choice + return true; + } else if (b > 0) { + // 'b' is greater than 0, so it is the better choice + return false; + } else { + // Both 'a' and 'b' are non-positive, prefer 'a' + return true; + } + }); + if (layerI != cumulativeEnergyInLayer.end()) { + int layer = std::distance(cumulativeEnergyInLayer.begin(), layerI); + if (TracksterLinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced) { + LogDebug("TracksterLinkingbySkeletons") + << "Layer containing at least an energy fraction " << percentage << " is " << layer << "\n"; + } + return rhtools_.getPositionLayer(layer, false).z(); + } else { + return 0.; + } +} + +void TracksterLinkingbySkeletons::linkTracksters( + const Inputs &input, + std::vector &resultTracksters, + std::vector> &linkedResultTracksters, + std::vector> &linkedTracksterIdToInputTracksterId) { + const auto &tracksters = input.tracksters; + const auto &layerClusters = input.layerClusters; + + // linking : trackster is hadr nic if its barycenter is in CE-H + auto isHadron = [&](const Trackster &t) -> bool { + auto boundary_z = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(); + + return (std::abs(t.barycenter().Z()) > boundary_z); + }; + + if (TracksterLinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced) + LogDebug("TracksterLinkingbySkeletons") << "------- Graph Linking ------- \n"; + + auto intersectLineWithSurface = [](float surfaceZ, const Vector &origin, const Vector &direction) -> Vector { + auto const t = (surfaceZ - origin.Z()) / direction.Z(); + auto const iX = t * direction.X() + origin.X(); + auto const iY = t * direction.Y() + origin.Y(); + auto const iZ = surfaceZ; + const Vector intersection(iX, iY, iZ); + return intersection; + }; + + auto returnSkeletons = [&](const Trackster &trackster) -> std::array { + auto const &vertices = trackster.vertices(); + std::vector vertices_lcs(vertices.size()); + std::transform(vertices.begin(), vertices.end(), vertices_lcs.begin(), [&layerClusters](unsigned int index) { + return layerClusters[index]; + }); + + std::sort(vertices_lcs.begin(), vertices_lcs.end(), [](reco::CaloCluster &c1, reco::CaloCluster &c2) { + return c1.position().z() < c2.position().z(); + }); + + auto const firstLayerZ = + findSkeletonPoints(0.1f, trackster.raw_energy(), trackster.vertices(), rhtools_, layerClusters); + auto const lastLayerZ = + findSkeletonPoints(0.9f, trackster.raw_energy(), trackster.vertices(), rhtools_, layerClusters); + auto const t0_p1 = trackster.barycenter(); + auto const t0_p0 = intersectLineWithSurface(firstLayerZ, t0_p1, trackster.eigenvectors(0)); + auto const t0_p2 = intersectLineWithSurface(lastLayerZ, t0_p1, trackster.eigenvectors(0)); + std::array skeleton{{t0_p0, t0_p1, t0_p2}}; + std::sort(skeleton.begin(), skeleton.end(), [](Vector &v1, Vector &v2) { return v1.Z() < v2.Z(); }); + return skeleton; + }; + + // auto argSortTrackster = [](const std::vector &tracksters) -> std::vector { + // std::vector retIndices(tracksters.size()); + // std::iota(retIndices.begin(), retIndices.end(), 0); + // std::stable_sort(retIndices.begin(), retIndices.end(), [&tracksters](unsigned int i, unsigned int j) { + // return tracksters[i].raw_energy() > tracksters[j].raw_energy(); + // }); + // + // return retIndices; + // }; + + auto isPointInCone = [](const Vector &origin, + const Vector &direction, + const float halfAngle, + const float maxHeight, + const Vector &testPoint, + std::string str = "") -> bool { + Vector toCheck = testPoint - origin; + auto projection = toCheck.Dot(direction.Unit()); + auto const angle = ROOT::Math::VectorUtil::Angle(direction, toCheck); + LogDebug("TracksterLinkingbySkeletons") + << str << "Origin " << origin << " TestPoint " << testPoint << " projection " << projection << " maxHeight " + << maxHeight << " Angle " << angle << " halfAngle " << halfAngle << std::endl; + if (projection < 0.f || projection > maxHeight) { + return false; + } + + return angle < halfAngle; + }; + + TICLLayerTile tracksterTilePos; + TICLLayerTile tracksterTileNeg; + + for (size_t id_t = 0; id_t < tracksters.size(); ++id_t) { + auto t = tracksters[id_t]; + if (t.barycenter().eta() > 0.) { + tracksterTilePos.fill(t.barycenter().eta(), t.barycenter().phi(), id_t); + } else if (t.barycenter().eta() < 0.) { + tracksterTileNeg.fill(t.barycenter().eta(), t.barycenter().phi(), id_t); + } + } + + //actual trackster-trackster linking + auto pcaQuality = [](const Trackster &trackster) -> float { + auto const &eigenvalues = trackster.eigenvalues(); + auto const e0 = eigenvalues[0]; + auto const e1 = eigenvalues[1]; + auto const e2 = eigenvalues[2]; + auto const sum = e0 + e1 + e2; + auto const normalized_e0 = e0 / sum; + auto const normalized_e1 = e1 / sum; + auto const normalized_e2 = e2 / sum; + return normalized_e0; + }; + + const float halfAngle0 = angle_first_cone_; + const float halfAngle1 = angle_second_cone_; + const float halfAngle2 = angle_third_cone_; + const float maxHeightCone = max_height_cone_; + + std::vector maskReceivedLink(tracksters.size(), 1); + std::vector isRootTracksters(tracksters.size(), 1); + + std::vector allNodes; + for (size_t it = 0; it < tracksters.size(); ++it) { + allNodes.emplace_back(it); + } + + for (size_t it = 0; it < tracksters.size(); ++it) { + auto const &trackster = tracksters[it]; + isHadron(trackster); + + auto pcaQ = pcaQuality(trackster); + LogDebug("TracksterLinkingbySkeletons") + << "DEBUG Trackster " << it << " energy " << trackster.raw_energy() << " Num verties " + << trackster.vertices().size() << " PCA Quality " << pcaQ << std::endl; + const float pcaQTh = pcaQ_; + const unsigned int pcaQLCSize = pcaQLCSize_; + if (pcaQ >= pcaQTh && trackster.vertices().size() > pcaQLCSize) { + auto const skeletons = returnSkeletons(trackster); + LogDebug("TracksterLinkingbySkeletons") + << "Trackster " << it << " energy " << trackster.raw_energy() << " Num verties " + << trackster.vertices().size() << " PCA Quality " << pcaQ << " Skeletons " << skeletons[0] << std::endl; + auto const &eigenVec = trackster.eigenvectors(0); + auto const eigenVal = trackster.eigenvalues()[0]; + auto const &directionOrigin = eigenVec * eigenVal; + + auto const bary = trackster.barycenter(); + float eta_min = std::max(abs(bary.eta()) - del_, TileConstants::minEta); + float eta_max = std::min(abs(bary.eta()) + del_, TileConstants::maxEta); + + if (bary.eta() > 0.) { + std::array search_box = + tracksterTilePos.searchBoxEtaPhi(eta_min, eta_max, bary.phi() - del_, bary.phi() + del_); + if (search_box[2] > search_box[3]) { + search_box[3] += TileConstants::nPhiBins; + } + + for (int eta_i = search_box[0]; eta_i <= search_box[1]; ++eta_i) { + for (int phi_i = search_box[2]; phi_i <= search_box[3]; ++phi_i) { + auto &neighbours = tracksterTilePos[tracksterTilePos.globalBin(eta_i, (phi_i % TileConstants::nPhiBins))]; + for (auto n : neighbours) { + if (maskReceivedLink[n] == 0) + continue; + + auto const &trackster_out = tracksters[n]; + auto const &skeletons_out = returnSkeletons(trackster_out); + auto const skeletonDist2 = (skeletons[2] - skeletons_out[0]).Mag2(); + auto const pcaQOuter = pcaQuality(trackster_out); + // auto const dotProd = ((skeletons[0]-skeletons[2]).Unit()).Dot((skeletons_out[0] - skeletons_out[2]).Unit()); + bool isGoodPCA = (pcaQOuter >= pcaQTh) && (trackster_out.vertices().size() > pcaQLCSize); + auto const maxHeightSmallCone = std::sqrt((skeletons[2] - skeletons[0]).Mag2()); + bool isInSmallCone = isPointInCone( + skeletons[0], directionOrigin, halfAngle0, maxHeightSmallCone, skeletons_out[0], "Small Cone"); + bool isInCone = + isPointInCone(skeletons[1], directionOrigin, halfAngle1, maxHeightCone, skeletons_out[0], "BigCone "); + bool isInLastCone = + isPointInCone(skeletons[2], directionOrigin, halfAngle2, maxHeightCone, skeletons_out[0], "LastCone"); + bool dotProd = + isGoodPCA + ? ((skeletons[0] - skeletons[2]).Unit()).Dot((skeletons_out[0] - skeletons_out[2]).Unit()) >= + dotCut_ + : true; + + LogDebug("TracksterLinkingbySkeletons") + << "\tTrying to Link Trackster " << n << " energy " << tracksters[n].raw_energy() << " LCs " + << tracksters[n].vertices().size() << " skeletons " << skeletons_out[0] << " Dist " << skeletonDist2 + << " dot Prod " + << ((skeletons[0] - skeletons[2]).Unit()).Dot((skeletons_out[0] - skeletons_out[2]).Unit()) + << " isGoodDotProd " << dotProd << " isPointInBigCone " << isInCone << " isPointInSmallCone " + << isInSmallCone << " isPointInLastCone " << isInLastCone << std::endl; + if (isInLastCone && dotProd) { + LogDebug("TracksterLinkingbySkeletons") << "\t==== LINK: Trackster " << it << " Linked with Trackster " + << n << " LCs " << tracksters[n].vertices().size() << std::endl; + LogDebug("TracksterLinkingbySkeletons") + << "\t\tSkeleton origin " << skeletons[2] << " Skeleton out " << skeletons_out[0] << std::endl; + maskReceivedLink[n] = 0; + allNodes[it].addNeighbour(n); + isRootTracksters[n] = 0; + } + if (isInCone && skeletonDist2 <= maxDistSkeletonsSq_ && dotProd) { + LogDebug("TracksterLinkingbySkeletons") << "\t==== LINK: Trackster " << it << " Linked with Trackster " + << n << " LCs " << tracksters[n].vertices().size() << std::endl; + LogDebug("TracksterLinkingbySkeletons") + << "\t\tSkeleton origin " << skeletons[2] << " Skeleton out " << skeletons_out[0] << std::endl; + maskReceivedLink[n] = 0; + allNodes[it].addNeighbour(n); + isRootTracksters[n] = 0; + continue; + } + if (isInSmallCone && !isGoodPCA) { + maskReceivedLink[n] = 0; + allNodes[it].addNeighbour(n); + isRootTracksters[n] = 0; + LogDebug("TracksterLinkingbySkeletons") + << "\t==== LINK: Trackster " << it << " Linked with Trackster in small cone " << n << std::endl; + LogDebug("TracksterLinkingbySkeletons") + << "\t\tSkeleton origin " << skeletons[0] << " Skeleton out " << skeletons_out[0] << std::endl; + continue; + } + } + } + } + } + } + } + + LogDebug("TracksterLinkingbySkeletons") << "**************** FINAL GRAPH **********************" << std::endl; + for (auto const &node : allNodes) { + if (isRootTracksters[node.getId()]) { + LogDebug("TracksterLinkingbySkeletons") + << "ISROOT " + << " Node " << node.getId() << " position " << tracksters[node.getId()].barycenter() << " energy " + << tracksters[node.getId()].raw_energy() << std::endl; + } else { + LogDebug("TracksterLinkingbySkeletons") + << "Node " << node.getId() << " position " << tracksters[node.getId()].barycenter() << " energy " + << tracksters[node.getId()].raw_energy() << std::endl; + } + } + LogDebug("TracksterLinkingbySkeletons") << "********************************************************" << std::endl; + + TICLGraph graph(allNodes, isRootTracksters); + + int ic = 0; + auto const &components = graph.findSubComponents(); + linkedTracksterIdToInputTracksterId.resize(components.size()); + for (auto const &comp : components) { + LogDebug("TracksterLinkingbySkeletons") << "Component " << ic << " Node: "; + std::vector linkedTracksters; + Trackster outTrackster; + for (auto const &node : comp) { + LogDebug("TracksterLinkingbySkeletons") << node << " "; + linkedTracksterIdToInputTracksterId[ic].push_back(node); + outTrackster.mergeTracksters(input.tracksters[node]); + } + linkedTracksters.push_back(resultTracksters.size()); + resultTracksters.push_back(outTrackster); + linkedResultTracksters.push_back(linkedTracksters); + LogDebug("TracksterLinkingbySkeletons") << "\n"; + ++ic; + } + LogDebug("TracksterLinkingbySkeletons") << "resultLinked " << std::endl; +} // linkTracksters diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.h b/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.h new file mode 100644 index 0000000000000..defbebb67856b --- /dev/null +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.h @@ -0,0 +1,89 @@ +#ifndef RecoHGCal_TICL_TracksterLinkingAlgoBySkeletons_H +#define RecoHGCal_TICL_TracksterLinkingAlgoBySkeletons_H + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "CommonTools/Utils/interface/StringCutObjectSelector.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" +#include "RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h" +#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "DataFormats/GeometrySurface/interface/BoundDisk.h" +#include "MagneticField/Engine/interface/MagneticField.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" + +namespace ticl { + + class TracksterLinkingbySkeletons : public TracksterLinkingAlgoBase { + public: + TracksterLinkingbySkeletons(const edm::ParameterSet& conf, edm::ConsumesCollector iC); + + virtual ~TracksterLinkingbySkeletons() {} + + void linkTracksters(const Inputs& input, + std::vector& resultTracksters, + std::vector>& linkedResultTracksters, + std::vector>& linkedTracksterIdToInputTracksterId) override; + + float findSkeletonPoints(float percentage, + const float trackster_energy, + const std::vector vertices, + const hgcal::RecHitTools& rhtools, + const std::vector& layerClusters); + + void initialize(const HGCalDDDConstants* hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) override; + + static void fillPSetDescription(edm::ParameterSetDescription& iDesc) { + iDesc.add("track_time_quality_threshold", 0.5); + iDesc.add("wind", 1.5); + iDesc.add("angle0", 1.523599); + iDesc.add("angle1", 1.349006); + iDesc.add("angle2", 1.174532); + iDesc.add("maxConeHeight", 500.); + iDesc.add("pcaQuality", 0.97); + iDesc.add("pcaQualityLCSize", 5); + iDesc.add("dotProdCut", 0.975); + iDesc.add("maxDistSkeletonsSq", 2500.); + TracksterLinkingAlgoBase::fillPSetDescription(iDesc); + } + + private: + using Vector = ticl::Trackster::Vector; + + void buildLayers(); + + void dumpLinksFound(std::vector>& resultCollection, const char* label) const; + + float timing_quality_threshold_; + float del_; + float angle_first_cone_; + float angle_second_cone_; + float angle_third_cone_; + float pcaQ_; + unsigned int pcaQLCSize_; + float dotCut_; + float maxDistSkeletonsSq_; + float max_height_cone_; + const HGCalDDDConstants* hgcons_; + + std::unique_ptr firstDisk_[2]; + std::unique_ptr interfaceDisk_[2]; + + hgcal::RecHitTools rhtools_; + + edm::ESHandle bfield_; + edm::ESHandle propagator_; + }; + +} // namespace ticl + +#endif diff --git a/RecoHGCal/TICL/plugins/TracksterLinksProducer.cc b/RecoHGCal/TICL/plugins/TracksterLinksProducer.cc index f9a53945a9e7f..72ef4d7623a05 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinksProducer.cc +++ b/RecoHGCal/TICL/plugins/TracksterLinksProducer.cc @@ -23,6 +23,17 @@ #include "RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h" #include "RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" + +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" + +#include "MagneticField/Engine/interface/MagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" + +#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" @@ -54,16 +65,29 @@ class TracksterLinksProducer : public edm::stream::EDProducer<> { std::vector>> original_masks_tokens_; - const edm::ESGetToken geometry_token_; std::once_flag initializeGeometry_; + + const edm::ESGetToken geometry_token_; + const std::string detector_; + const std::string propName_; + + const edm::ESGetToken bfield_token_; + const edm::ESGetToken propagator_token_; + const HGCalDDDConstants *hgcons_; hgcal::RecHitTools rhtools_; + edm::ESGetToken hdc_token_; }; TracksterLinksProducer::TracksterLinksProducer(const edm::ParameterSet &ps) : clusters_token_(consumes>(ps.getParameter("layer_clusters"))), clustersTime_token_( consumes>>(ps.getParameter("layer_clustersTime"))), - geometry_token_(esConsumes()) { + geometry_token_(esConsumes()), + detector_(ps.getParameter("detector")), + propName_(ps.getParameter("propagator")), + bfield_token_(esConsumes()), + propagator_token_( + esConsumes(edm::ESInputTag("", propName_))) { // Loop over the edm::VInputTag and append the token to tracksters_tokens_ for (auto const &tag : ps.getParameter>("tracksters_collections")) { tracksters_tokens_.emplace_back(consumes>(tag)); @@ -84,6 +108,13 @@ TracksterLinksProducer::TracksterLinksProducer(const edm::ParameterSet &ps) auto linkingPSet = ps.getParameter("linkingPSet"); auto algoType = linkingPSet.getParameter("type"); + + if (algoType == "Skeletons") { + std::string detectorName_ = (detector_ == "HFNose") ? "HGCalHFNoseSensitive" : "HGCalEESensitive"; + hdc_token_ = esConsumes( + edm::ESInputTag("", detectorName_)); + } + linkingAlgo_ = TracksterLinkingPluginFactory::get()->create(algoType, linkingPSet, consumesCollector()); } @@ -92,8 +123,16 @@ void TracksterLinksProducer::beginJob() {} void TracksterLinksProducer::endJob(){}; void TracksterLinksProducer::beginRun(edm::Run const &iEvent, edm::EventSetup const &es) { + edm::ESHandle hdc = es.getHandle(hdc_token_); + hgcons_ = hdc.product(); + edm::ESHandle geom = es.getHandle(geometry_token_); rhtools_.setGeometry(*geom); + + edm::ESHandle bfield = es.getHandle(bfield_token_); + edm::ESHandle propagator = es.getHandle(propagator_token_); + + linkingAlgo_->initialize(hgcons_, rhtools_, bfield, propagator); }; void TracksterLinksProducer::dumpTrackster(const Trackster &t) const { @@ -203,7 +242,8 @@ void TracksterLinksProducer::printTrackstersDebug(const std::vector & void TracksterLinksProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { edm::ParameterSetDescription desc; edm::ParameterSetDescription linkingDesc; - linkingDesc.addNode(edm::PluginDescription("type", "FastJet", true)); + linkingDesc.addNode(edm::PluginDescription("type", "Skeletons", true)); + desc.add("linkingPSet", linkingDesc); desc.add>( "tracksters_collections", {edm::InputTag("ticlTrackstersCLUE3DEM"), edm::InputTag("ticlTrackstersCLUE3DHAD")}); @@ -211,6 +251,8 @@ void TracksterLinksProducer::fillDescriptions(edm::ConfigurationDescriptions &de {edm::InputTag("hgcalMergeLayerClusters", "InitialLayerClustersMask")}); desc.add("layer_clusters", edm::InputTag("hgcalMergeLayerClusters")); desc.add("layer_clustersTime", edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster")); + desc.add("detector", "HGCAL"); + desc.add("propagator", "PropagatorWithMaterial"); descriptions.add("tracksterLinksProducer", desc); } diff --git a/RecoHGCal/TICL/python/customiseForTICLv5_cff.py b/RecoHGCal/TICL/python/customiseForTICLv5_cff.py index 60f1979e5e251..4f55f544ee641 100644 --- a/RecoHGCal/TICL/python/customiseForTICLv5_cff.py +++ b/RecoHGCal/TICL/python/customiseForTICLv5_cff.py @@ -14,6 +14,9 @@ from RecoHGCal.TICL.ticlCandidateProducer_cfi import ticlCandidateProducer as _ticlCandidateProducer from RecoHGCal.Configuration.RecoHGCal_EventContent_cff import customiseForTICLv5EventContent from RecoHGCal.TICL.iterativeTICL_cff import ticlIterLabels, ticlIterLabelsMerge +from SimCalorimetry.HGCalAssociatorProducers.TSToSimTSAssociation_cfi import tracksterSimTracksterAssociationLinkingbyCLUE3D as _tracksterSimTracksterAssociationLinkingbyCLUE3D +from SimCalorimetry.HGCalAssociatorProducers.TSToSimTSAssociation_cfi import tracksterSimTracksterAssociationPRbyCLUE3D as _tracksterSimTracksterAssociationPRbyCLUE3D +from Validation.HGCalValidation.HGCalValidator_cff import hgcalValidatorv5 def customiseForTICLv5(process): @@ -31,12 +34,52 @@ def customiseForTICLv5(process): process.ticlCandidate = _ticlCandidateProducer.clone() process.ticlCandidateTask = cms.Task(process.ticlCandidate) + process.tracksterSimTracksterAssociationLinkingbyCLUE3DEM = _tracksterSimTracksterAssociationLinkingbyCLUE3D.clone( + label_tst = cms.InputTag("ticlTrackstersCLUE3DEM") + ) + process.tracksterSimTracksterAssociationPRbyCLUE3DEM = _tracksterSimTracksterAssociationPRbyCLUE3D.clone( + label_tst = cms.InputTag("ticlTrackstersCLUE3DEM") + ) + process.tracksterSimTracksterAssociationLinkingbyCLUE3DHAD = _tracksterSimTracksterAssociationLinkingbyCLUE3D.clone( + label_tst = cms.InputTag("ticlTrackstersCLUE3DHAD") + ) + process.tracksterSimTracksterAssociationPRbyCLUE3DHAD = _tracksterSimTracksterAssociationPRbyCLUE3D.clone( + label_tst = cms.InputTag("ticlTrackstersCLUE3DHAD") + ) + process.iterTICLTask = cms.Task(process.ticlLayerTileTask, process.ticlIterationsTask, process.ticlTracksterLinksTask, process.ticlCandidateTask) + process.particleFlowClusterHGCal.initialClusteringStep.tracksterSrc = "ticlCandidate" + process.globalrecoTask.remove(process.ticlTrackstersMerge) + + process.tracksterSimTracksterAssociationLinking.label_tst = cms.InputTag("ticlCandidate") + process.tracksterSimTracksterAssociationPR.label_tst = cms.InputTag("ticlCandidate") + + process.tracksterSimTracksterAssociationLinkingPU.label_tst = cms.InputTag("ticlTracksterLinks") + process.tracksterSimTracksterAssociationPRPU.label_tst = cms.InputTag("ticlTracksterLinks") process.mergeTICLTask = cms.Task() + process.pfTICL.ticlCandidateSrc = cms.InputTag("ticlCandidate") + process.hgcalAssociators = cms.Task(process.lcAssocByEnergyScoreProducer, process.layerClusterCaloParticleAssociationProducer, + process.scAssocByEnergyScoreProducer, process.layerClusterSimClusterAssociationProducer, + process.lcSimTSAssocByEnergyScoreProducer, process.layerClusterSimTracksterAssociationProducer, + process.simTsAssocByEnergyScoreProducer, process.simTracksterHitLCAssociatorByEnergyScoreProducer, + process.tracksterSimTracksterAssociationLinking, process.tracksterSimTracksterAssociationPR, + process.tracksterSimTracksterAssociationLinkingbyCLUE3DEM, process.tracksterSimTracksterAssociationPRbyCLUE3DEM, + process.tracksterSimTracksterAssociationLinkingbyCLUE3DHAD, process.tracksterSimTracksterAssociationPRbyCLUE3DHAD, + process.tracksterSimTracksterAssociationLinkingPU, process.tracksterSimTracksterAssociationPRPU + ) + + process.hgcalValidatorv5 = hgcalValidatorv5.clone( + ticlTrackstersMerge = cms.InputTag("ticlCandidate"), + trackstersclue3d = cms.InputTag("ticlTracksterLinks") + ) + process.hgcalValidatorSequence = cms.Sequence(process.hgcalValidatorv5) + process.hgcalValidation = cms.Sequence(process.hgcalSimHitValidationEE+process.hgcalSimHitValidationHEF+process.hgcalSimHitValidationHEB+process.hgcalDigiValidationEE+process.hgcalDigiValidationHEF+process.hgcalDigiValidationHEB+process.hgcalRecHitValidationEE+process.hgcalRecHitValidationHEF+process.hgcalRecHitValidationHEB+process.hgcalHitValidationSequence+process.hgcalValidatorSequence+process.hgcalTiclPFValidation+process.hgcalPFJetValidation) + process.globalValidationHGCal = cms.Sequence(process.hgcalValidation) + process.validation_step9 = cms.EndPath(process.globalValidationHGCal) process = customiseForTICLv5EventContent(process) diff --git a/Validation/HGCalValidation/python/HGCalValidator_cfi.py b/Validation/HGCalValidation/python/HGCalValidator_cfi.py index eae93bf63017e..391818e0b5615 100644 --- a/Validation/HGCalValidation/python/HGCalValidator_cfi.py +++ b/Validation/HGCalValidation/python/HGCalValidator_cfi.py @@ -14,6 +14,7 @@ labelTst.extend([cms.InputTag("ticlSimTracksters", "fromCPs"), cms.InputTag("ticlSimTracksters")]) lcInputMask = [cms.InputTag("ticlTracksters"+iteration) for iteration in ticlIterLabels] lcInputMask.extend([cms.InputTag("ticlSimTracksters", "fromCPs"), cms.InputTag("ticlSimTracksters")]) + hgcalValidator = DQMEDAnalyzer( "HGCalValidator", @@ -91,6 +92,90 @@ ) +#labelTst = [cms.InputTag("ticlTracksters"+iteration) for iteration in ticlIterLabelsMerge] +labelTst = ["ticlTrackstersCLUE3DEM", "ticlTrackstersCLUE3DHAD", "ticlTracksterLinks"] +labelTst.extend([cms.InputTag("ticlSimTracksters", "fromCPs"), cms.InputTag("ticlSimTracksters")]) +lcInputMask = ["ticlTrackstersCLUE3DEM", "ticlTrackstersCLUE3DHAD", "ticlTracksterLinks"] +#lcInputMask = [cms.InputTag("ticlTracksters"+iteration) for iteration in ticlIterLabels] +lcInputMask.extend([cms.InputTag("ticlSimTracksters", "fromCPs"), cms.InputTag("ticlSimTracksters")]) + +hgcalValidatorv5 = DQMEDAnalyzer( + "HGCalValidator", + + ### general settings ### + # selection of CP for evaluation of efficiency # + CaloParticleSelectionForEfficiency, + + ### reco input configuration ### + #2DLayerClusters, PFClusters, Tracksters + label_lcl = layerClusterCaloParticleAssociation.label_lc, + label_tst = cms.VInputTag(labelTst), + label_simTS = cms.InputTag("ticlSimTracksters"), + label_simTSFromCP = cms.InputTag("ticlSimTracksters", "fromCPs"), + + associator = cms.untracked.InputTag("layerClusterCaloParticleAssociationProducer"), + + associatorSim = cms.untracked.InputTag("layerClusterSimClusterAssociationProducer"), + + #General info on layers etc. + SaveGeneralInfo = cms.untracked.bool(True), + #CaloParticle related plots + doCaloParticlePlots = cms.untracked.bool(True), + #Select caloParticles for efficiency or pass through + doCaloParticleSelection = cms.untracked.bool(True), + #SimCluster related plots + doSimClustersPlots = cms.untracked.bool(True), + label_SimClusters = cms.InputTag("SimClusters"), + label_SimClustersLevel = cms.InputTag("ClusterLevel"), + #Layer Cluster related plots + doLayerClustersPlots = cms.untracked.bool(True), + label_layerClusterPlots = cms.InputTag("hgcalMergeLayerClusters"), + label_LCToCPLinking = cms.InputTag("LCToCP_association"), + #Trackster related plots + doTrackstersPlots = cms.untracked.bool(True), + label_TS = cms.string("Morphology"), + label_TSToCPLinking = cms.string("TSToCP_linking"), + label_TSToSTSPR = cms.string("TSToSTS_patternRecognition"), + #candidates plots + doCandidatesPlots = cms.untracked.bool(True), + ticlCandidates = cms.InputTag("ticlCandidate"), + + ticlTrackstersMerge = cms.InputTag("ticlTrackstersMerge"), + simTiclCandidates = cms.InputTag("ticlSimTracksters"), + recoTracks = cms.InputTag("generalTracks"), + trackstersclue3d = cms.InputTag("ticlCandidate"), + mergeRecoToSimAssociator = cms.InputTag("tracksterSimTracksterAssociationLinking", "recoToSim"), + mergeSimToRecoAssociator = cms.InputTag("tracksterSimTracksterAssociationLinking", "simToReco"), + + #The cumulative material budget in front of each layer. To be more specific, it + #is the material budget just in front of the active material (not including it). + #This file is created using the official material budget code. + cummatbudinxo = cms.FileInPath('Validation/HGCalValidation/data/D41.cumulative.xo'), + + ### sim input configuration ### + label_cp_effic = layerClusterCaloParticleAssociation.label_cp, + label_cp_fake = cms.InputTag("mix","MergedCaloTruth"), + #simClusters + label_scl = layerClusterSimClusterAssociation.label_scl, + + simVertices = cms.InputTag("g4SimHits"), + + LayerClustersInputMask = cms.VInputTag(lcInputMask), + + #Total number of layers of HGCal that we want to monitor + #Could get this also from HGCalImagingAlgo::maxlayer but better to get it from here + totallayers_to_monitor = cms.int32(52), + #Thicknesses we want to monitor. -1 is for scintillator + thicknesses_to_monitor = cms.vint32(120,200,300,-1), + + # HistoProducerAlgo. Defines the set of plots to be booked and filled + histoProducerAlgoBlock = HGVHistoProducerAlgoBlock, + + ### output configuration + dirName = cms.string('HGCAL/HGCalValidator/') + +) + from Configuration.ProcessModifiers.premix_stage2_cff import premix_stage2 premix_stage2.toModify(hgcalValidator, label_cp_fake = "mixData:MergedCaloTruth"