CMS 3D CMS Logo

List of all members | Classes | Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes
ticl::PatternRecognitionbyCLUE3D< TILES > Class Template Referencefinal

#include <PatternRecognitionbyCLUE3D.h>

Inheritance diagram for ticl::PatternRecognitionbyCLUE3D< TILES >:
ticl::PatternRecognitionAlgoBaseT< TILES >

Classes

struct  ClustersOnLayer
 

Public Member Functions

void energyRegressionAndID (const std::vector< reco::CaloCluster > &layerClusters, const tensorflow::Session *, std::vector< Trackster > &result)
 
void makeTracksters (const typename PatternRecognitionAlgoBaseT< TILES >::Inputs &input, std::vector< Trackster > &result, std::unordered_map< int, std::vector< int >> &seedToTracksterAssociation) override
 
 PatternRecognitionbyCLUE3D (const edm::ParameterSet &conf, edm::ConsumesCollector)
 
 ~PatternRecognitionbyCLUE3D () override=default
 
- Public Member Functions inherited from ticl::PatternRecognitionAlgoBaseT< TILES >
virtual void makeTracksters (const Inputs &input, std::vector< Trackster > &result, std::unordered_map< int, std::vector< int >> &seedToTracksterAssociation)=0
 
 PatternRecognitionAlgoBaseT (const edm::ParameterSet &conf, edm::ConsumesCollector)
 
virtual ~PatternRecognitionAlgoBaseT ()
 

Static Public Member Functions

static void fillPSetDescription (edm::ParameterSetDescription &iDesc)
 

Private Member Functions

void calculateDistanceToHigher (const TILES &, const int layerId, const std::vector< std::pair< int, int >> &)
 
void calculateLocalDensity (const TILES &, const int layerId, const std::vector< std::pair< int, int >> &)
 
void dumpClusters (const TILES &tiles, const std::vector< std::pair< int, int >> &layerIdx2layerandSoa, const int) const
 
void dumpTiles (const TILES &) const
 
void dumpTracksters (const std::vector< std::pair< int, int >> &layerIdx2layerandSoa, const int, const std::vector< Trackster > &) const
 
int findAndAssignTracksters (const TILES &, const std::vector< std::pair< int, int >> &)
 
void reset ()
 

Private Attributes

edm::ESGetToken< CaloGeometry, CaloGeometryRecordcaloGeomToken_
 
std::vector< ClustersOnLayerclusters_
 
const bool computeLocalTime_
 
const std::vector< double > criticalDensity_
 
const std::vector< double > criticalEtaPhiDistance_
 
const std::vector< double > criticalSelfDensity_
 
const std::vector< double > criticalXYDistance_
 
const std::vector< int > criticalZDistanceLyr_
 
const float cutHadProb_
 
const std::vector< double > densityEtaPhiDistanceSqr_
 
const bool densityOnSameLayer_
 
const std::vector< int > densitySiblingLayers_
 
const std::vector< double > densityXYDistanceSqr_
 
const bool doPidCut_
 
const std::string eidInputName_
 
const float eidMinClusterEnergy_
 
const int eidNClusters_
 
const int eidNLayers_
 
const std::string eidOutputNameEnergy_
 
const std::string eidOutputNameId_
 
tensorflow::Session * eidSession_
 
const std::vector< int > filter_on_categories_
 
const std::vector< double > kernelDensityFactor_
 
std::vector< float > layersPosZ_
 
const std::vector< int > minNumLayerCluster_
 
const bool nearestHigherOnSameLayer_
 
const std::vector< double > outlierMultiplier_
 
const bool rescaleDensityByZ_
 
hgcal::RecHitTools rhtools_
 
std::vector< int > tracksterSeedAlgoId_
 
const bool useAbsoluteProjectiveScale_
 
const bool useClusterDimensionXY_
 
const bool usePCACleaning_
 

Static Private Attributes

static const int eidNFeatures_ = 3
 

Additional Inherited Members

- Protected Attributes inherited from ticl::PatternRecognitionAlgoBaseT< TILES >
int algo_verbosity_
 

Detailed Description

template<typename TILES>
class ticl::PatternRecognitionbyCLUE3D< TILES >

Definition at line 12 of file PatternRecognitionbyCLUE3D.h.

Constructor & Destructor Documentation

◆ PatternRecognitionbyCLUE3D()

template<typename TILES >
PatternRecognitionbyCLUE3D::PatternRecognitionbyCLUE3D ( const edm::ParameterSet conf,
edm::ConsumesCollector  iC 
)

Definition at line 21 of file PatternRecognitionbyCLUE3D.cc.

24  criticalDensity_(conf.getParameter<std::vector<double>>("criticalDensity")),
25  criticalSelfDensity_(conf.getParameter<std::vector<double>>("criticalSelfDensity")),
26  densitySiblingLayers_(conf.getParameter<std::vector<int>>("densitySiblingLayers")),
27  densityEtaPhiDistanceSqr_(conf.getParameter<std::vector<double>>("densityEtaPhiDistanceSqr")),
28  densityXYDistanceSqr_(conf.getParameter<std::vector<double>>("densityXYDistanceSqr")),
29  kernelDensityFactor_(conf.getParameter<std::vector<double>>("kernelDensityFactor")),
30  densityOnSameLayer_(conf.getParameter<bool>("densityOnSameLayer")),
31  nearestHigherOnSameLayer_(conf.getParameter<bool>("nearestHigherOnSameLayer")),
32  useAbsoluteProjectiveScale_(conf.getParameter<bool>("useAbsoluteProjectiveScale")),
33  useClusterDimensionXY_(conf.getParameter<bool>("useClusterDimensionXY")),
34  rescaleDensityByZ_(conf.getParameter<bool>("rescaleDensityByZ")),
35  criticalEtaPhiDistance_(conf.getParameter<std::vector<double>>("criticalEtaPhiDistance")),
36  criticalXYDistance_(conf.getParameter<std::vector<double>>("criticalXYDistance")),
37  criticalZDistanceLyr_(conf.getParameter<std::vector<int>>("criticalZDistanceLyr")),
38  outlierMultiplier_(conf.getParameter<std::vector<double>>("outlierMultiplier")),
39  minNumLayerCluster_(conf.getParameter<std::vector<int>>("minNumLayerCluster")),
40  doPidCut_(conf.getParameter<bool>("doPidCut")),
41  cutHadProb_(conf.getParameter<double>("cutHadProb")),
42  eidInputName_(conf.getParameter<std::string>("eid_input_name")),
43  eidOutputNameEnergy_(conf.getParameter<std::string>("eid_output_name_energy")),
44  eidOutputNameId_(conf.getParameter<std::string>("eid_output_name_id")),
45  eidMinClusterEnergy_(conf.getParameter<double>("eid_min_cluster_energy")),
46  eidNLayers_(conf.getParameter<int>("eid_n_layers")),
47  eidNClusters_(conf.getParameter<int>("eid_n_clusters")),
48  computeLocalTime_(conf.getParameter<bool>("computeLocalTime")),
49  usePCACleaning_(conf.getParameter<bool>("usePCACleaning")){};
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeomToken_
const std::vector< double > densityEtaPhiDistanceSqr_
const std::vector< double > densityXYDistanceSqr_
const std::vector< double > outlierMultiplier_
const std::vector< int > densitySiblingLayers_
const std::vector< double > kernelDensityFactor_
const std::vector< double > criticalDensity_
const std::vector< double > criticalXYDistance_
const std::vector< double > criticalEtaPhiDistance_
const std::vector< int > criticalZDistanceLyr_
const std::vector< double > criticalSelfDensity_

◆ ~PatternRecognitionbyCLUE3D()

template<typename TILES >
ticl::PatternRecognitionbyCLUE3D< TILES >::~PatternRecognitionbyCLUE3D ( )
overridedefault

Member Function Documentation

◆ calculateDistanceToHigher()

template<typename TILES >
void PatternRecognitionbyCLUE3D::calculateDistanceToHigher ( const TILES &  tiles,
const int  layerId,
const std::vector< std::pair< int, int >> &  layerIdx2layerandSoa 
)
private

Definition at line 684 of file PatternRecognitionbyCLUE3D.cc.

References funct::abs(), ticl::Advanced, hltPFPuppi_cfi::algoId, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), HLT_2024v14_cff::delta_phi, reco::deltaPhi(), reco::deltaR2(), mps_fire::i, hcalRecHitTable_cff::ieta, hcalRecHitTable_cff::iphi, SiStripPI::max, allConversions_cfi::maxDelta, SiStripPI::min, HLT_2024v14_cff::minLayer, L1TMuonDQMOffline_cfi::nEtaBins, ecaldqm::binning::nPhiBins, electrons_cff::numberOfClusters, hltrates_dqm_sourceclient-live_cfg::offset, trackingPOGFilters_cfi::phiWindow, and mathSSE::sqrt().

685  {
688  auto &clustersOnLayer = clusters_[layerId];
689  unsigned int numberOfClusters = clustersOnLayer.x.size();
690 
691  auto distanceSqr = [](float r0, float r1, float phi0, float phi1) -> float {
692  auto delta_phi = reco::deltaPhi(phi0, phi1);
693  return (r0 - r1) * (r0 - r1) + r1 * r1 * delta_phi * delta_phi;
694  };
695 
696  for (unsigned int i = 0; i < numberOfClusters; i++) {
698  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
699  << "Starting searching nearestHigher on " << layerId << " with rho: " << clustersOnLayer.rho[i]
700  << " at eta, phi: " << tiles[layerId].etaBin(clustersOnLayer.eta[i]) << ", "
701  << tiles[layerId].phiBin(clustersOnLayer.phi[i]);
702  }
703  // We need to partition the two sides of the HGCAL detector
704  auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
705  int minLayer = 0;
706  int maxLayer = 2 * lastLayerPerSide - 1;
707  auto algoId = clustersOnLayer.algoId[i];
708  if (layerId < lastLayerPerSide) {
710  maxLayer = std::min(layerId + densitySiblingLayers_[algoId], lastLayerPerSide - 1);
711  } else {
712  minLayer = std::max(layerId - densitySiblingLayers_[algoId], lastLayerPerSide + 1);
713  maxLayer = std::min(layerId + densitySiblingLayers_[algoId], maxLayer);
714  }
716  float i_delta = maxDelta;
717  std::pair<int, int> i_nearestHigher(-1, -1);
718  std::pair<float, int> nearest_distances(maxDelta, std::numeric_limits<int>::max());
719  for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) {
720  if (!nearestHigherOnSameLayer_ && (layerId == currentLayer))
721  continue;
722  const auto &tileOnLayer = tiles[currentLayer];
723  int etaWindow = 1;
724  int phiWindow = 1;
725  int etaBinMin = std::max(tileOnLayer.etaBin(clustersOnLayer.eta[i]) - etaWindow, 0);
726  int etaBinMax = std::min(tileOnLayer.etaBin(clustersOnLayer.eta[i]) + etaWindow, nEtaBin - 1);
727  int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow;
728  int phiBinMax = tileOnLayer.phiBin(clustersOnLayer.phi[i]) + phiWindow;
729  for (int ieta = etaBinMin; ieta <= etaBinMax; ++ieta) {
730  auto offset = ieta * nPhiBin;
731  for (int iphi_it = phiBinMin; iphi_it <= phiBinMax; ++iphi_it) {
732  int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin);
734  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
735  << "Searching nearestHigher on " << currentLayer << " eta, phi: " << ieta << ", " << iphi_it << " "
736  << iphi << " " << offset << " " << (offset + iphi);
737  }
738  for (auto otherClusterIdx : tileOnLayer[offset + iphi]) {
739  auto const &layerandSoa = layerIdx2layerandSoa[otherClusterIdx];
740  // Skip masked layer clusters
741  if ((layerandSoa.first == -1) && (layerandSoa.second == -1))
742  continue;
743  auto const &clustersOnOtherLayer = clusters_[layerandSoa.first];
744  auto dist = maxDelta;
745  auto dist_transverse = maxDelta;
746  int dist_layers = std::abs(layerandSoa.first - layerId);
748  dist_transverse = distanceSqr(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
749  clustersOnOtherLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
750  clustersOnLayer.phi[i],
751  clustersOnOtherLayer.phi[layerandSoa.second]);
752  // Add Z-scale to the final distance
753  dist = dist_transverse;
754  } else {
755  dist = reco::deltaR2(clustersOnLayer.eta[i],
756  clustersOnLayer.phi[i],
757  clustersOnOtherLayer.eta[layerandSoa.second],
758  clustersOnOtherLayer.phi[layerandSoa.second]);
759  dist_transverse = dist;
760  }
761  bool foundHigher = (clustersOnOtherLayer.rho[layerandSoa.second] > clustersOnLayer.rho[i]) ||
762  (clustersOnOtherLayer.rho[layerandSoa.second] == clustersOnLayer.rho[i] &&
763  clustersOnOtherLayer.layerClusterOriginalIdx[layerandSoa.second] >
764  clustersOnLayer.layerClusterOriginalIdx[i]);
766  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
767  << "Searching nearestHigher on " << currentLayer
768  << " with rho: " << clustersOnOtherLayer.rho[layerandSoa.second]
769  << " on layerIdxInSOA: " << layerandSoa.first << ", " << layerandSoa.second
770  << " with distance: " << sqrt(dist) << " foundHigher: " << foundHigher;
771  }
772  if (foundHigher && dist <= i_delta) {
773  // update i_delta
774  i_delta = dist;
775  nearest_distances = std::make_pair(sqrt(dist_transverse), dist_layers);
776  // update i_nearestHigher
777  i_nearestHigher = layerandSoa;
778  }
779  } // End of loop on clusters
780  } // End of loop on phi bins
781  } // End of loop on eta bins
782  } // End of loop on layers
783 
784  clustersOnLayer.delta[i] = nearest_distances;
785  clustersOnLayer.nearestHigher[i] = i_nearestHigher;
786  } // End of loop on clusters
787 }
constexpr double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:26
Log< level::Info, true > LogVerbatim
std::vector< ClustersOnLayer > clusters_
const std::vector< int > densitySiblingLayers_
T sqrt(T t)
Definition: SSEVec.h:23
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
static constexpr int nPhiBins
unsigned int lastLayer(bool nose=false) const
Definition: RecHitTools.h:80

◆ calculateLocalDensity()

template<typename TILES >
void PatternRecognitionbyCLUE3D::calculateLocalDensity ( const TILES &  tiles,
const int  layerId,
const std::vector< std::pair< int, int >> &  layerIdx2layerandSoa 
)
private

Definition at line 518 of file PatternRecognitionbyCLUE3D.cc.

References funct::abs(), ticl::Advanced, hltPFPuppi_cfi::algoId, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), HLT_2024v14_cff::delta_phi, reco::deltaPhi(), reco::deltaR2(), f, mps_fire::i, hcalRecHitTable_cff::ieta, ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets::if(), hcalRecHitTable_cff::iphi, SiStripPI::max, SiStripPI::min, HLT_2024v14_cff::minLayer, L1TMuonDQMOffline_cfi::nEtaBins, ecaldqm::binning::nPhiBins, electrons_cff::numberOfClusters, hltrates_dqm_sourceclient-live_cfg::offset, trackingPOGFilters_cfi::phiWindow, mathSSE::sqrt(), testProducerWithPsetDescEmpty_cfi::y1, and testProducerWithPsetDescEmpty_cfi::y2.

519  {
522  auto &clustersOnLayer = clusters_[layerId];
523  unsigned int numberOfClusters = clustersOnLayer.x.size();
524 
525  auto isReachable = [](float r0, float r1, float phi0, float phi1, float delta_sqr) -> bool {
526  auto delta_phi = reco::deltaPhi(phi0, phi1);
527  return (r0 - r1) * (r0 - r1) + r1 * r1 * delta_phi * delta_phi < delta_sqr;
528  };
529  auto distance_debug = [&](float x1, float x2, float y1, float y2) -> float {
530  return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
531  };
532 
533  for (unsigned int i = 0; i < numberOfClusters; i++) {
534  auto algoId = clustersOnLayer.algoId[i];
535  // We need to partition the two sides of the HGCAL detector
536  auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
537  int minLayer = 0;
538  int maxLayer = 2 * lastLayerPerSide - 1;
539  if (layerId < lastLayerPerSide) {
541  maxLayer = std::min(layerId + densitySiblingLayers_[algoId], lastLayerPerSide - 1);
542  } else {
543  minLayer = std::max(layerId - densitySiblingLayers_[algoId], lastLayerPerSide);
544  maxLayer = std::min(layerId + densitySiblingLayers_[algoId], maxLayer);
545  }
546  float deltaLayersZ = std::abs(layersPosZ_[maxLayer % lastLayerPerSide] - layersPosZ_[minLayer % lastLayerPerSide]);
547 
548  for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) {
550  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "RefLayer: " << layerId << " SoaIDX: " << i;
551  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "NextLayer: " << currentLayer;
552  }
553  const auto &tileOnLayer = tiles[currentLayer];
554  bool onSameLayer = (currentLayer == layerId);
556  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "onSameLayer: " << onSameLayer;
557  }
558  const int etaWindow = 2;
559  const int phiWindow = 2;
560  int etaBinMin = std::max(tileOnLayer.etaBin(clustersOnLayer.eta[i]) - etaWindow, 0);
561  int etaBinMax = std::min(tileOnLayer.etaBin(clustersOnLayer.eta[i]) + etaWindow, nEtaBin - 1);
562  int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow;
563  int phiBinMax = tileOnLayer.phiBin(clustersOnLayer.phi[i]) + phiWindow;
565  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "eta: " << clustersOnLayer.eta[i];
566  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "phi: " << clustersOnLayer.phi[i];
567  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "etaBinMin: " << etaBinMin << ", etaBinMax: " << etaBinMax;
568  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "phiBinMin: " << phiBinMin << ", phiBinMax: " << phiBinMax;
569  }
570  for (int ieta = etaBinMin; ieta <= etaBinMax; ++ieta) {
571  auto offset = ieta * nPhiBin;
573  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "offset: " << offset;
574  }
575  for (int iphi_it = phiBinMin; iphi_it <= phiBinMax; ++iphi_it) {
576  int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin);
578  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "iphi: " << iphi;
579  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
580  << "Entries in tileBin: " << tileOnLayer[offset + iphi].size();
581  }
582  for (auto otherClusterIdx : tileOnLayer[offset + iphi]) {
583  auto const &layerandSoa = layerIdx2layerandSoa[otherClusterIdx];
584  // Skip masked layer clusters
585  if ((layerandSoa.first == -1) && (layerandSoa.second == -1)) {
587  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Skipping masked layerIdx " << otherClusterIdx;
588  }
589  continue;
590  }
591  auto const &clustersLayer = clusters_[layerandSoa.first];
593  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
594  << "OtherLayer: " << layerandSoa.first << " SoaIDX: " << layerandSoa.second;
595  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherEta: " << clustersLayer.eta[layerandSoa.second];
596  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherPhi: " << clustersLayer.phi[layerandSoa.second];
597  }
598 
599  bool onSameCluster = clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx;
600  if (onSameLayer && !densityOnSameLayer_ && !onSameCluster) {
602  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
603  << "Skipping different cluster " << otherClusterIdx << "in the same layer " << currentLayer;
604  }
605  continue;
606  }
607 
608  bool reachable = false;
611  reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
612  clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
613  clustersOnLayer.phi[i],
614  clustersLayer.phi[layerandSoa.second],
615  clustersOnLayer.radius[i] * clustersOnLayer.radius[i]);
616  } else {
617  // Still differentiate between silicon and Scintillator.
618  // Scintillator has yet to be studied further.
619  if (clustersOnLayer.isSilicon[i]) {
620  reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
621  clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
622  clustersOnLayer.phi[i],
623  clustersLayer.phi[layerandSoa.second],
625  } else {
626  reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
627  clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
628  clustersOnLayer.phi[i],
629  clustersLayer.phi[layerandSoa.second],
630  clustersOnLayer.radius[i] * clustersOnLayer.radius[i]);
631  }
632  }
633  } else {
634  reachable = (reco::deltaR2(clustersOnLayer.eta[i],
635  clustersOnLayer.phi[i],
636  clustersLayer.eta[layerandSoa.second],
637  clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_[algoId]);
638  }
640  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[eta,phi]: "
641  << reco::deltaR2(clustersOnLayer.eta[i],
642  clustersOnLayer.phi[i],
643  clustersLayer.eta[layerandSoa.second],
644  clustersLayer.phi[layerandSoa.second]);
645  auto dist = distance_debug(
646  clustersOnLayer.r_over_absz[i],
647  clustersLayer.r_over_absz[layerandSoa.second],
648  clustersOnLayer.r_over_absz[i] * std::abs(clustersOnLayer.phi[i]),
649  clustersLayer.r_over_absz[layerandSoa.second] * std::abs(clustersLayer.phi[layerandSoa.second]));
650  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[cm]: " << (dist * clustersOnLayer.z[i]);
651  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
652  << "Energy Other: " << clustersLayer.energy[layerandSoa.second];
653  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Cluster radius: " << clustersOnLayer.radius[i];
654  }
655  if (reachable) {
656  float factor_same_layer_different_cluster = (onSameLayer && !densityOnSameLayer_) ? 0.f : 1.f;
657  auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx
658  ? 1.f
659  : kernelDensityFactor_[algoId] * factor_same_layer_different_cluster) *
660  clustersLayer.energy[layerandSoa.second];
661  clustersOnLayer.rho[i] += energyToAdd;
662  clustersOnLayer.z_extension[i] = deltaLayersZ;
664  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
665  << "Adding " << energyToAdd << " partial " << clustersOnLayer.rho[i];
666  }
667  }
668  } // end of loop on possible compatible clusters
669  } // end of loop over phi-bin region
670  } // end of loop over eta-bin region
671  } // end of loop on the sibling layers
672  if (rescaleDensityByZ_) {
674  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
675  << "Rescaling original density: " << clustersOnLayer.rho[i] << " by Z: " << deltaLayersZ
676  << " to final density/cm: " << clustersOnLayer.rho[i] / deltaLayersZ;
677  }
678  clustersOnLayer.rho[i] /= deltaLayersZ;
679  }
680  } // end of loop over clusters on this layer
681 }
constexpr double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:26
Log< level::Info, true > LogVerbatim
const std::vector< double > densityEtaPhiDistanceSqr_
const std::vector< double > densityXYDistanceSqr_
std::vector< ClustersOnLayer > clusters_
const std::vector< int > densitySiblingLayers_
T sqrt(T t)
Definition: SSEVec.h:23
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const std::vector< double > kernelDensityFactor_
double f[11][100]
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
static constexpr int nPhiBins
if(threadIdxLocalY==0 &&threadIdxLocalX==0)
unsigned int lastLayer(bool nose=false) const
Definition: RecHitTools.h:80

◆ dumpClusters()

template<typename TILES >
void PatternRecognitionbyCLUE3D::dumpClusters ( const TILES &  tiles,
const std::vector< std::pair< int, int >> &  layerIdx2layerandSoa,
const int  eventNumber 
) const
private

Definition at line 109 of file PatternRecognitionbyCLUE3D.cc.

References ticl::Advanced, alignBH_cfg::fixed, EgammaValidation_cff::num, and findQualityFiles::v.

111  {
113  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "[evt, lyr, Seed, x, y, z, r/|z|, eta, phi, "
114  "etab, phib, cells, enrgy, e/rho, rho, z_ext, "
115  " dlt_tr, dlt_lyr, "
116  " nestHL, nestHSoaIdx, radius, clIdx, lClOrigIdx, SOAidx";
117  }
118 
119  for (unsigned int layer = 0; layer < clusters_.size(); layer++) {
120  auto const &thisLayer = clusters_[layer];
121  int num = 0;
122  for (auto v : thisLayer.x) {
124  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
125  << std::setw(4) << eventNumber << ", " << std::setw(3) << layer << ", " << std::setw(4)
126  << thisLayer.isSeed[num] << ", " << std::setprecision(3) << std::fixed << v << ", " << thisLayer.y[num]
127  << ", " << thisLayer.z[num] << ", " << thisLayer.r_over_absz[num] << ", " << thisLayer.eta[num] << ", "
128  << thisLayer.phi[num] << ", " << std::setw(5) << tiles[layer].etaBin(thisLayer.eta[num]) << ", "
129  << std::setw(5) << tiles[layer].phiBin(thisLayer.phi[num]) << ", " << std::setw(4) << thisLayer.cells[num]
130  << ", " << std::setprecision(3) << thisLayer.energy[num] << ", "
131  << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " << thisLayer.rho[num] << ", "
132  << thisLayer.z_extension[num] << ", " << std::scientific << thisLayer.delta[num].first << ", "
133  << std::setw(10) << thisLayer.delta[num].second << ", " << std::setw(5)
134  << thisLayer.nearestHigher[num].first << ", " << std::setw(10) << thisLayer.nearestHigher[num].second
135  << ", " << std::defaultfloat << std::setprecision(3) << thisLayer.radius[num] << ", " << std::setw(5)
136  << thisLayer.clusterIndex[num] << ", " << std::setw(4) << thisLayer.layerClusterOriginalIdx[num] << ", "
137  << std::setw(4) << num << ", ClusterInfo";
138  }
139  ++num;
140  }
141  }
142  for (unsigned int lcIdx = 0; lcIdx < layerIdx2layerandSoa.size(); lcIdx++) {
143  auto const &layerandSoa = layerIdx2layerandSoa[lcIdx];
144  // Skip masked layer clusters
145  if ((layerandSoa.first == -1) && (layerandSoa.second == -1))
146  continue;
148  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
149  << "lcIdx: " << lcIdx << " on Layer: " << layerandSoa.first << " SOA: " << layerandSoa.second;
150  }
151  }
152 }
Log< level::Info, true > LogVerbatim
std::vector< ClustersOnLayer > clusters_

◆ dumpTiles()

template<typename TILES >
void PatternRecognitionbyCLUE3D::dumpTiles ( const TILES &  tiles) const
private

Definition at line 52 of file PatternRecognitionbyCLUE3D.cc.

References ticl::Advanced, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), relativeConstraints::empty, hcalRecHitTable_cff::ieta, hcalRecHitTable_cff::iphi, L1TMuonDQMOffline_cfi::nEtaBins, ecaldqm::binning::nPhiBins, and hltrates_dqm_sourceclient-live_cfg::offset.

52  {
55  auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
56  int maxLayer = 2 * lastLayerPerSide - 1;
57  for (int layer = 0; layer <= maxLayer; layer++) {
58  for (int ieta = 0; ieta < nEtaBin; ieta++) {
59  auto offset = ieta * nPhiBin;
60  for (int phi = 0; phi < nPhiBin; phi++) {
61  int iphi = ((phi % nPhiBin + nPhiBin) % nPhiBin);
62  if (!tiles[layer][offset + iphi].empty()) {
64  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer: " << layer << " ieta: " << ieta << " phi: " << phi
65  << " " << tiles[layer][offset + iphi].size();
66  }
67  }
68  }
69  }
70  }
71 }
Log< level::Info, true > LogVerbatim
static constexpr int nPhiBins
unsigned int lastLayer(bool nose=false) const
Definition: RecHitTools.h:80

◆ dumpTracksters()

template<typename TILES >
void PatternRecognitionbyCLUE3D::dumpTracksters ( const std::vector< std::pair< int, int >> &  layerIdx2layerandSoa,
const int  eventNumber,
const std::vector< Trackster > &  tracksters 
) const
private

Definition at line 74 of file PatternRecognitionbyCLUE3D.cc.

References ticl::Advanced, ticl::Trackster::charged_hadron, ticl::Trackster::electron, ticl::Trackster::neutral_hadron, EgammaValidation_cff::num, ticl::Trackster::photon, AlCaHLTBitMon_QueryRunRegistry::string, submitPVValidationJobs::t, and findQualityFiles::v.

76  {
78  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
79  << "[evt, tracksterId, cells, prob_photon, prob_ele, prob_chad, prob_nhad, layer_i, x_i, y_i, eta_i, phi_i, "
80  "energy_i, radius_i, rho_i, z_extension, delta_tr, delta_lyr, isSeed_i";
81  }
82 
83  int num = 0;
84  const std::string sep(", ");
85  for (auto const &t : tracksters) {
86  for (auto v : t.vertices()) {
87  auto [lyrIdx, soaIdx] = layerIdx2layerandSoa[v];
88  auto const &thisLayer = clusters_[lyrIdx];
90  edm::LogVerbatim("PatternRecognitionbyCLUE3D_NTP")
91  << std::setw(4) << eventNumber << sep << std::setw(4) << num << sep << std::setw(4) << t.vertices().size()
92  << sep << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::photon) << sep << std::setw(8)
93  << t.id_probability(ticl::Trackster::ParticleType::electron) << sep << std::setw(8)
94  << t.id_probability(ticl::Trackster::ParticleType::charged_hadron) << sep << std::setw(8)
95  << t.id_probability(ticl::Trackster::ParticleType::neutral_hadron) << sep << std::setw(4) << lyrIdx << sep
96  << std::setw(10) << thisLayer.x[soaIdx] << sep << std::setw(10) << thisLayer.y[soaIdx] << sep
97  << std::setw(10) << thisLayer.eta[soaIdx] << sep << std::setw(10) << thisLayer.phi[soaIdx] << sep
98  << std::setw(10) << thisLayer.energy[soaIdx] << sep << std::setw(10) << thisLayer.radius[soaIdx] << sep
99  << std::setw(10) << thisLayer.rho[soaIdx] << sep << std::setw(10) << thisLayer.z_extension[soaIdx] << sep
100  << std::setw(10) << thisLayer.delta[soaIdx].first << sep << std::setw(10) << thisLayer.delta[soaIdx].second
101  << sep << std::setw(4) << thisLayer.isSeed[soaIdx];
102  }
103  }
104  num++;
105  }
106 }
Log< level::Info, true > LogVerbatim
std::vector< ClustersOnLayer > clusters_

◆ energyRegressionAndID()

template<typename TILES >
void PatternRecognitionbyCLUE3D::energyRegressionAndID ( const std::vector< reco::CaloCluster > &  layerClusters,
const tensorflow::Session *  eidSession,
std::vector< Trackster > &  result 
)

Definition at line 378 of file PatternRecognitionbyCLUE3D.cc.

References a, funct::abs(), b, HLT_FULL_cff::batchSize, data, hcalRecHitTable_cff::energy, reco::CaloCluster::energy(), reco::CaloCluster::eta(), f, lowptgsfeleseed::features(), nano_mu_digi_cff::float, reco::CaloCluster::hitsAndFractions(), mps_fire::i, input, crabTemplate::inputList, createfilelist::int, dqmiolumiharvest::j, dqmdumpme::k, MainPageGenerator::l, hltEgammaHGCALIDVarsL1Seeded_cfi::layerClusters, jetsAK4_CHS_cff::outputNames, PatBasicFWLiteJetAnalyzer_Selector_cfg::outputs, reco::CaloCluster::phi(), tensorflow::run(), l1trig_cff::shape, jetUpdater_cfi::sort, bphysicsOniaDQM_cfi::vertex, ticl::Trackster::vertex_multiplicity(), AlignmentTracksFromVertexSelector_cfi::vertices, and ticl::Trackster::vertices().

380  {
381  // Energy regression and particle identification strategy:
382  //
383  // 1. Set default values for regressed energy and particle id for each trackster.
384  // 2. Store indices of tracksters whose total sum of cluster energies is above the
385  // eidMinClusterEnergy_ (GeV) threshold. Inference is not applied for soft tracksters.
386  // 3. When no trackster passes the selection, return.
387  // 4. Create input and output tensors. The batch dimension is determined by the number of
388  // selected tracksters.
389  // 5. Fill input tensors with layer cluster features. Per layer, clusters are ordered descending
390  // by energy. Given that tensor data is contiguous in memory, we can use pointer arithmetic to
391  // fill values, even with batching.
392  // 6. Zero-fill features for empty clusters in each layer.
393  // 7. Batched inference.
394  // 8. Assign the regressed energy and id probabilities to each trackster.
395  //
396  // Indices used throughout this method:
397  // i -> batch element / trackster
398  // j -> layer
399  // k -> cluster
400  // l -> feature
401 
402  // set default values per trackster, determine if the cluster energy threshold is passed,
403  // and store indices of hard tracksters
404  std::vector<int> tracksterIndices;
405  for (int i = 0; i < static_cast<int>(tracksters.size()); i++) {
406  // calculate the cluster energy sum (2)
407  // note: after the loop, sumClusterEnergy might be just above the threshold which is enough to
408  // decide whether to run inference for the trackster or not
409  float sumClusterEnergy = 0.;
410  for (const unsigned int &vertex : tracksters[i].vertices()) {
411  sumClusterEnergy += static_cast<float>(layerClusters[vertex].energy());
412  // there might be many clusters, so try to stop early
413  if (sumClusterEnergy >= eidMinClusterEnergy_) {
414  // set default values (1)
415  tracksters[i].setRegressedEnergy(0.f);
416  tracksters[i].zeroProbabilities();
417  tracksterIndices.push_back(i);
418  break;
419  }
420  }
421  }
422 
423  // do nothing when no trackster passes the selection (3)
424  int batchSize = static_cast<int>(tracksterIndices.size());
425  if (batchSize == 0) {
426  return;
427  }
428 
429  // create input and output tensors (4)
430  tensorflow::TensorShape shape({batchSize, eidNLayers_, eidNClusters_, eidNFeatures_});
431  tensorflow::Tensor input(tensorflow::DT_FLOAT, shape);
433 
434  std::vector<tensorflow::Tensor> outputs;
435  std::vector<std::string> outputNames;
436  if (!eidOutputNameEnergy_.empty()) {
438  }
439  if (!eidOutputNameId_.empty()) {
440  outputNames.push_back(eidOutputNameId_);
441  }
442 
443  // fill input tensor (5)
444  for (int i = 0; i < batchSize; i++) {
445  const Trackster &trackster = tracksters[tracksterIndices[i]];
446 
447  // per layer, we only consider the first eidNClusters_ clusters in terms of energy, so in order
448  // to avoid creating large / nested structures to do the sorting for an unknown number of total
449  // clusters, create a sorted list of layer cluster indices to keep track of the filled clusters
450  std::vector<int> clusterIndices(trackster.vertices().size());
451  for (int k = 0; k < (int)trackster.vertices().size(); k++) {
452  clusterIndices[k] = k;
453  }
454  sort(clusterIndices.begin(), clusterIndices.end(), [&layerClusters, &trackster](const int &a, const int &b) {
455  return layerClusters[trackster.vertices(a)].energy() > layerClusters[trackster.vertices(b)].energy();
456  });
457 
458  // keep track of the number of seen clusters per layer
459  std::vector<int> seenClusters(eidNLayers_);
460 
461  // loop through clusters by descending energy
462  for (const int &k : clusterIndices) {
463  // get features per layer and cluster and store the values directly in the input tensor
464  const reco::CaloCluster &cluster = layerClusters[trackster.vertices(k)];
465  int j = rhtools_.getLayerWithOffset(cluster.hitsAndFractions()[0].first) - 1;
466  if (j < eidNLayers_ && seenClusters[j] < eidNClusters_) {
467  // get the pointer to the first feature value for the current batch, layer and cluster
468  float *features = &input.tensor<float, 4>()(i, j, seenClusters[j], 0);
469 
470  // fill features
471  *(features++) = float(cluster.energy() / float(trackster.vertex_multiplicity(k)));
472  *(features++) = float(std::abs(cluster.eta()));
473  *(features) = float(cluster.phi());
474 
475  // increment seen clusters
476  seenClusters[j]++;
477  }
478  }
479 
480  // zero-fill features of empty clusters in each layer (6)
481  for (int j = 0; j < eidNLayers_; j++) {
482  for (int k = seenClusters[j]; k < eidNClusters_; k++) {
483  float *features = &input.tensor<float, 4>()(i, j, k, 0);
484  for (int l = 0; l < eidNFeatures_; l++) {
485  *(features++) = 0.f;
486  }
487  }
488  }
489  }
490 
491  // run the inference (7)
493 
494  // store regressed energy per trackster (8)
495  if (!eidOutputNameEnergy_.empty()) {
496  // get the pointer to the energy tensor, dimension is batch x 1
497  float *energy = outputs[0].flat<float>().data();
498 
499  for (const int &i : tracksterIndices) {
500  tracksters[i].setRegressedEnergy(*(energy++));
501  }
502  }
503 
504  // store id probabilities per trackster (8)
505  if (!eidOutputNameId_.empty()) {
506  // get the pointer to the id probability tensor, dimension is batch x id_probabilities.size()
507  int probsIdx = eidOutputNameEnergy_.empty() ? 0 : 1;
508  float *probs = outputs[probsIdx].flat<float>().data();
509 
510  for (const int &i : tracksterIndices) {
511  tracksters[i].setProbabilities(probs);
512  probs += tracksters[i].id_probabilities().size();
513  }
514  }
515 }
std::vector< NamedTensor > NamedTensorList
Definition: TensorFlow.h:31
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:210
double phi() const
azimuthal angle of cluster centroid
Definition: CaloCluster.h:184
static std::string const input
Definition: EdmProvDump.cc:50
std::vector< float > features(const reco::PreId &ecal, const reco::PreId &hcal, double rho, const reco::BeamSpot &spot, noZS::EcalClusterLazyTools &ecalTools)
void run(Session *session, const NamedTensorList &inputs, const std::vector< std::string > &outputNames, std::vector< Tensor > *outputs, const thread::ThreadPoolOptions &threadPoolOptions)
Definition: TensorFlow.cc:271
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double f[11][100]
string inputList
Definition: crabTemplate.py:6
double energy() const
cluster energy
Definition: CaloCluster.h:149
std::vector< unsigned int > & vertices()
Definition: Trackster.h:57
std::vector< float > & vertex_multiplicity()
Definition: Trackster.h:58
double b
Definition: hdecay.h:120
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:80
double a
Definition: hdecay.h:121
double eta() const
pseudorapidity of cluster centroid
Definition: CaloCluster.h:181
unsigned int getLayerWithOffset(const DetId &) const
Definition: RecHitTools.cc:381

◆ fillPSetDescription()

template<typename TILES >
void PatternRecognitionbyCLUE3D::fillPSetDescription ( edm::ParameterSetDescription iDesc)
static

Definition at line 863 of file PatternRecognitionbyCLUE3D.cc.

References edm::ParameterSetDescription::add(), and AlCaHLTBitMon_QueryRunRegistry::string.

863  {
864  iDesc.add<int>("algo_verbosity", 0);
865  iDesc.add<std::vector<double>>("criticalDensity", {4, 4, 4})->setComment("in GeV");
866  iDesc.add<std::vector<double>>("criticalSelfDensity", {0.15, 0.15, 0.15} /* roughly 1/(densitySiblingLayers+1) */)
867  ->setComment("Minimum ratio of self_energy/local_density to become a seed.");
868  iDesc.add<std::vector<int>>("densitySiblingLayers", {3, 3, 3})
869  ->setComment(
870  "inclusive, layers to consider while computing local density and searching for nearestHigher higher");
871  iDesc.add<std::vector<double>>("densityEtaPhiDistanceSqr", {0.0008, 0.0008, 0.0008})
872  ->setComment("in eta,phi space, distance to consider for local density");
873  iDesc.add<std::vector<double>>("densityXYDistanceSqr", {3.24, 3.24, 3.24})
874  ->setComment("in cm, distance on the transverse plane to consider for local density");
875  iDesc.add<std::vector<double>>("kernelDensityFactor", {0.2, 0.2, 0.2})
876  ->setComment("Kernel factor to be applied to other LC while computing the local density");
877  iDesc.add<bool>("densityOnSameLayer", false);
878  iDesc.add<bool>("nearestHigherOnSameLayer", false)
879  ->setComment("Allow the nearestHigher to be located on the same layer");
880  iDesc.add<bool>("useAbsoluteProjectiveScale", true)
881  ->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables");
882  iDesc.add<bool>("useClusterDimensionXY", false)
883  ->setComment(
884  "Boolean. If true use the estimated cluster radius to determine the cluster compatibility while computing "
885  "the local density");
886  iDesc.add<bool>("rescaleDensityByZ", false)
887  ->setComment(
888  "Rescale local density by the extension of the Z 'volume' explored. The transvere dimension is, at present, "
889  "fixed and factored out.");
890  iDesc.add<std::vector<double>>("criticalEtaPhiDistance", {0.025, 0.025, 0.025})
891  ->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed");
892  iDesc.add<std::vector<double>>("criticalXYDistance", {1.8, 1.8, 1.8})
893  ->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed");
894  iDesc.add<std::vector<int>>("criticalZDistanceLyr", {5, 5, 5})
895  ->setComment("Minimal distance in layers along the Z axis from nearestHigher to become a seed");
896  iDesc.add<std::vector<double>>("outlierMultiplier", {2, 2, 2})
897  ->setComment("Minimal distance in transverse space from nearestHigher to become an outlier");
898  iDesc.add<std::vector<int>>("minNumLayerCluster", {2, 2, 2})->setComment("Not Inclusive");
899  iDesc.add<bool>("doPidCut", false);
900  iDesc.add<double>("cutHadProb", 0.5);
901  iDesc.add<std::string>("eid_input_name", "input");
902  iDesc.add<std::string>("eid_output_name_energy", "output/regressed_energy");
903  iDesc.add<std::string>("eid_output_name_id", "output/id_probabilities");
904  iDesc.add<double>("eid_min_cluster_energy", 1.);
905  iDesc.add<int>("eid_n_layers", 50);
906  iDesc.add<int>("eid_n_clusters", 10);
907  iDesc.add<bool>("computeLocalTime", false);
908  iDesc.add<bool>("usePCACleaning", false)->setComment("Enable PCA cleaning alorithm");
909 }
ParameterDescriptionBase * add(U const &iLabel, T const &value)

◆ findAndAssignTracksters()

template<typename TILES >
int PatternRecognitionbyCLUE3D::findAndAssignTracksters ( const TILES &  tiles,
const std::vector< std::pair< int, int >> &  layerIdx2layerandSoa 
)
private

Definition at line 790 of file PatternRecognitionbyCLUE3D.cc.

References ticl::Advanced, hltPFPuppi_cfi::algoId, mps_fire::i, and electrons_cff::numberOfClusters.

791  {
792  unsigned int nTracksters = 0;
793  std::vector<std::pair<int, int>> localStack;
794  const auto &critical_transverse_distance =
796  // find cluster seeds and outlier
797  for (unsigned int layer = 0; layer < 2 * rhtools_.lastLayer(); layer++) {
798  auto &clustersOnLayer = clusters_[layer];
799  unsigned int numberOfClusters = clustersOnLayer.x.size();
800  for (unsigned int i = 0; i < numberOfClusters; i++) {
801  // initialize clusterIndex
802  clustersOnLayer.clusterIndex[i] = -1;
803  auto algoId = clustersOnLayer.algoId[i];
804  bool isSeed = (clustersOnLayer.delta[i].first > critical_transverse_distance[algoId] ||
805  clustersOnLayer.delta[i].second > criticalZDistanceLyr_[algoId]) &&
806  (clustersOnLayer.rho[i] >= criticalDensity_[algoId]) &&
807  (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_[algoId]);
808  if (!clustersOnLayer.isSilicon[i]) {
809  isSeed = (clustersOnLayer.delta[i].first > clustersOnLayer.radius[i] ||
810  clustersOnLayer.delta[i].second > criticalZDistanceLyr_[algoId]) &&
811  (clustersOnLayer.rho[i] >= criticalDensity_[algoId]) &&
812  (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_[algoId]);
813  }
814  bool isOutlier =
815  (clustersOnLayer.delta[i].first > outlierMultiplier_[algoId] * critical_transverse_distance[algoId]) &&
816  (clustersOnLayer.rho[i] < criticalDensity_[algoId]);
817  if (isSeed) {
819  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
820  << "Found seed on Layer " << layer << " SOAidx: " << i << " assigned ClusterIdx: " << nTracksters;
821  }
822  clustersOnLayer.clusterIndex[i] = nTracksters++;
823  tracksterSeedAlgoId_.push_back(algoId);
824  clustersOnLayer.isSeed[i] = true;
825  localStack.emplace_back(layer, i);
826  } else if (!isOutlier) {
827  auto [lyrIdx, soaIdx] = clustersOnLayer.nearestHigher[i];
829  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
830  << "Found follower on Layer " << layer << " SOAidx: " << i << " attached to cluster on layer: " << lyrIdx
831  << " SOAidx: " << soaIdx;
832  }
833  if (lyrIdx >= 0)
834  clusters_[lyrIdx].followers[soaIdx].emplace_back(layer, i);
835  } else {
837  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
838  << "Found Outlier on Layer " << layer << " SOAidx: " << i << " with rho: " << clustersOnLayer.rho[i]
839  << " and delta: " << clustersOnLayer.delta[i].first << ", " << clustersOnLayer.delta[i].second;
840  }
841  }
842  }
843  }
844 
845  // Propagate cluster index
846  while (!localStack.empty()) {
847  auto [lyrIdx, soaIdx] = localStack.back();
848  auto &thisSeed = clusters_[lyrIdx].followers[soaIdx];
849  localStack.pop_back();
850 
851  // loop over followers
852  for (auto [follower_lyrIdx, follower_soaIdx] : thisSeed) {
853  // pass id to a follower
854  clusters_[follower_lyrIdx].clusterIndex[follower_soaIdx] = clusters_[lyrIdx].clusterIndex[soaIdx];
855  // push this follower to localStack
856  localStack.emplace_back(follower_lyrIdx, follower_soaIdx);
857  }
858  }
859  return nTracksters;
860 }
Log< level::Info, true > LogVerbatim
std::vector< ClustersOnLayer > clusters_
const std::vector< double > outlierMultiplier_
const std::vector< double > criticalDensity_
const std::vector< double > criticalXYDistance_
const std::vector< double > criticalEtaPhiDistance_
const std::vector< int > criticalZDistanceLyr_
unsigned int lastLayer(bool nose=false) const
Definition: RecHitTools.h:80
const std::vector< double > criticalSelfDensity_

◆ makeTracksters()

template<typename TILES >
void PatternRecognitionbyCLUE3D::makeTracksters ( const typename PatternRecognitionAlgoBaseT< TILES >::Inputs input,
std::vector< Trackster > &  result,
std::unordered_map< int, std::vector< int >> &  seedToTracksterAssociation 
)
override

Definition at line 155 of file PatternRecognitionbyCLUE3D.cc.

References funct::abs(), ticl::Advanced, cms::cuda::assert(), ticl::assignPCAtoTracksters(), ticl::Trackster::charged_hadron, hcalRecHitTable_cff::detId, f, relativeConstraints::geom, edm::EventSetup::getData(), reco::CaloCluster::hgcal_em, mps_fire::i, input, createfilelist::int, SiStripPI::max, ticl::Trackster::neutral_hadron, point, reset(), mps_fire::result, mathSSE::sqrt(), makeMuonMisalignmentScenario::sum_x, makeMuonMisalignmentScenario::sum_y, submitPVValidationJobs::t, findQualityFiles::v, and x.

158  {
159  // Protect from events with no seeding regions
160  if (input.regions.empty())
161  return;
162 
163  const int eventNumber = input.ev.eventAuxiliary().event();
165  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "New Event";
166  }
167 
168  edm::EventSetup const &es = input.es;
171 
172  // Assume identical Z-positioning between positive and negative sides.
173  // Also, layers inside the HGCAL geometry start from 1.
174  for (unsigned int i = 0; i < rhtools_.lastLayer(); ++i) {
175  layersPosZ_.push_back(rhtools_.getPositionLayer(i + 1).z());
177  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer " << i << " located at Z: " << layersPosZ_.back();
178  }
179  }
180 
181  clusters_.clear();
182  tracksterSeedAlgoId_.clear();
183 
184  clusters_.resize(2 * rhtools_.lastLayer(false));
185  std::vector<std::pair<int, int>> layerIdx2layerandSoa; //used everywhere also to propagate cluster masking
186 
187  layerIdx2layerandSoa.reserve(input.layerClusters.size());
188  unsigned int layerIdx = 0;
189  for (auto const &lc : input.layerClusters) {
190  if (input.mask[layerIdx] == 0.) {
192  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Skipping masked cluster: " << layerIdx;
193  }
194  layerIdx2layerandSoa.emplace_back(-1, -1);
195  layerIdx++;
196  continue;
197  }
198  const auto firstHitDetId = lc.hitsAndFractions()[0].first;
199  int layer = rhtools_.getLayerWithOffset(firstHitDetId) - 1 +
200  rhtools_.lastLayer(false) * ((rhtools_.zside(firstHitDetId) + 1) >> 1);
201  assert(layer >= 0);
202  auto detId = lc.hitsAndFractions()[0].first;
203  int layerClusterIndexInLayer = clusters_[layer].x.size();
204  layerIdx2layerandSoa.emplace_back(layer, layerClusterIndexInLayer);
205  float sum_x = 0.;
206  float sum_y = 0.;
207  float sum_sqr_x = 0.;
208  float sum_sqr_y = 0.;
209  float ref_x = lc.x();
210  float ref_y = lc.y();
211  float invClsize = 1. / lc.hitsAndFractions().size();
212  for (auto const &hitsAndFractions : lc.hitsAndFractions()) {
213  auto const &point = rhtools_.getPosition(hitsAndFractions.first);
214  sum_x += point.x() - ref_x;
215  sum_sqr_x += (point.x() - ref_x) * (point.x() - ref_x);
216  sum_y += point.y() - ref_y;
217  sum_sqr_y += (point.y() - ref_y) * (point.y() - ref_y);
218  }
219  // The variance of X for X uniform in circle of radius R is R^2/4,
220  // therefore we multiply the sqrt(var) by 2 to have a rough estimate of the
221  // radius. On the other hand, while averaging the x and y radius, we would
222  // end up dividing by 2. Hence we omit the value here and in the average
223  // below, too.
224  float radius_x = sqrt((sum_sqr_x - (sum_x * sum_x) * invClsize) * invClsize);
225  float radius_y = sqrt((sum_sqr_y - (sum_y * sum_y) * invClsize) * invClsize);
227  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
228  << "cluster rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5) << radius_y
229  << ", r: " << std::setw(5) << (radius_x + radius_y) << ", cells: " << std::setw(4)
230  << lc.hitsAndFractions().size();
231  }
232 
233  // The case of single cell layer clusters has to be handled differently.
234 
235  if (invClsize == 1.) {
236  // Silicon case
237  if (rhtools_.isSilicon(detId)) {
238  radius_x = radius_y = rhtools_.getRadiusToSide(detId);
240  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster in silicon, rx: " << std::setw(5)
241  << radius_x << ", ry: " << std::setw(5) << radius_y;
242  }
243  } else {
244  auto const &point = rhtools_.getPosition(detId);
245  auto const &eta_phi_window = rhtools_.getScintDEtaDPhi(detId);
246  radius_x = radius_y = point.perp() * eta_phi_window.second;
248  edm::LogVerbatim("PatternRecognitionbyCLUE3D")
249  << "Single cell cluster in scintillator. rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5)
250  << radius_y << ", eta-span: " << std::setw(5) << eta_phi_window.first << ", phi-span: " << std::setw(5)
251  << eta_phi_window.second;
252  }
253  }
254  }
255 
256  // Maybe check if these vectors can be reserved beforehands
257  clusters_[layer].x.emplace_back(lc.x());
258  clusters_[layer].y.emplace_back(lc.y());
259  clusters_[layer].z.emplace_back(lc.z());
260  clusters_[layer].r_over_absz.emplace_back(sqrt(lc.x() * lc.x() + lc.y() * lc.y()) / std::abs(lc.z()));
261  clusters_[layer].radius.emplace_back(radius_x + radius_y);
262  clusters_[layer].eta.emplace_back(lc.eta());
263  clusters_[layer].phi.emplace_back(lc.phi());
264  clusters_[layer].cells.push_back(lc.hitsAndFractions().size());
265  clusters_[layer].algoId.push_back(lc.algo() - reco::CaloCluster::hgcal_em);
266  clusters_[layer].isSilicon.push_back(rhtools_.isSilicon(detId));
267  clusters_[layer].energy.emplace_back(lc.energy());
268  clusters_[layer].isSeed.push_back(false);
269  clusters_[layer].clusterIndex.emplace_back(-1);
270  clusters_[layer].layerClusterOriginalIdx.emplace_back(layerIdx++);
271  clusters_[layer].nearestHigher.emplace_back(-1, -1);
272  clusters_[layer].rho.emplace_back(0.f);
273  clusters_[layer].z_extension.emplace_back(0.f);
274  clusters_[layer].delta.emplace_back(
276  }
277  for (unsigned int layer = 0; layer < clusters_.size(); layer++) {
278  clusters_[layer].followers.resize(clusters_[layer].x.size());
279  }
280 
281  auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
282  int maxLayer = 2 * lastLayerPerSide - 1;
283  std::vector<int> numberOfClustersPerLayer(maxLayer, 0);
284  for (int i = 0; i <= maxLayer; i++) {
285  calculateLocalDensity(input.tiles, i, layerIdx2layerandSoa);
286  }
287  for (int i = 0; i <= maxLayer; i++) {
288  calculateDistanceToHigher(input.tiles, i, layerIdx2layerandSoa);
289  }
290 
291  auto nTracksters = findAndAssignTracksters(input.tiles, layerIdx2layerandSoa);
293  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Reconstructed " << nTracksters << " tracksters" << std::endl;
294  dumpClusters(input.tiles, layerIdx2layerandSoa, eventNumber);
295  }
296 
297  // Build Trackster
298  result.resize(nTracksters);
299 
300  for (unsigned int layer = 0; layer < clusters_.size(); ++layer) {
301  const auto &thisLayer = clusters_[layer];
303  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Examining Layer: " << layer;
304  }
305  for (unsigned int lc = 0; lc < thisLayer.x.size(); ++lc) {
307  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Trackster " << thisLayer.clusterIndex[lc];
308  }
309  if (thisLayer.clusterIndex[lc] >= 0) {
311  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << " adding lcIdx: " << thisLayer.layerClusterOriginalIdx[lc];
312  }
313  result[thisLayer.clusterIndex[lc]].vertices().push_back(thisLayer.layerClusterOriginalIdx[lc]);
314  result[thisLayer.clusterIndex[lc]].vertex_multiplicity().push_back(1);
315  // loop over followers
316  for (auto [follower_lyrIdx, follower_soaIdx] : thisLayer.followers[lc]) {
317  std::array<unsigned int, 2> edge = {
318  {(unsigned int)thisLayer.layerClusterOriginalIdx[lc],
319  (unsigned int)clusters_[follower_lyrIdx].layerClusterOriginalIdx[follower_soaIdx]}};
320  result[thisLayer.clusterIndex[lc]].edges().push_back(edge);
321  }
322  }
323  }
324  }
325  size_t tracksterIndex = 0;
326  result.erase(std::remove_if(std::begin(result),
327  std::end(result),
328  [&](auto const &v) {
329  return static_cast<int>(v.vertices().size()) <
330  minNumLayerCluster_[tracksterSeedAlgoId_[tracksterIndex++]];
331  }),
332  result.end());
333  if (doPidCut_) {
334  energyRegressionAndID(input.layerClusters, input.tfSession, result);
335  result.erase(std::remove_if(std::begin(result),
336  std::end(result),
337  [&](auto const &v) {
338  auto const &hadProb =
341  return hadProb >= cutHadProb_;
342  }),
343  result.end());
344  }
345  result.shrink_to_fit();
346 
348  input.layerClusters,
349  input.layerClustersTime,
351  rhtools_,
353  true, // energy weighting
355 
357  for (auto const &t : result) {
358  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Barycenter: " << t.barycenter();
359  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "LCs: " << t.vertices().size();
360  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Energy: " << t.raw_energy();
361  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Regressed: " << t.regressed_energy();
362  }
363  }
364 
365  // Dump Tracksters information
367  dumpTracksters(layerIdx2layerandSoa, eventNumber, result);
368  }
369 
370  // Reset internal clusters_ structure of array for next event
371  reset();
373  edm::LogVerbatim("PatternRecognitionbyCLUE3D") << std::endl;
374  }
375 }
Log< level::Info, true > LogVerbatim
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeomToken_
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
void energyRegressionAndID(const std::vector< reco::CaloCluster > &layerClusters, const tensorflow::Session *, std::vector< Trackster > &result)
T z() const
Definition: PV3DBase.h:61
std::float_t getRadiusToSide(const DetId &) const
Definition: RecHitTools.cc:246
void assignPCAtoTracksters(std::vector< Trackster > &tracksters, const std::vector< reco::CaloCluster > &layerClusters, const edm::ValueMap< std::pair< float, float >> &layerClustersTime, double z_limit_em, hgcal::RecHitTools const &rhTools, bool computeLocalTime=false, bool energyWeight=true, bool clean=false, int minLayer=10, int maxLayer=10)
std::vector< ClustersOnLayer > clusters_
std::pair< float, float > getScintDEtaDPhi(const DetId &) const
Definition: RecHitTools.cc:236
assert(be >=bs)
int zside(const DetId &id) const
Definition: RecHitTools.cc:174
void calculateDistanceToHigher(const TILES &, const int layerId, const std::vector< std::pair< int, int >> &)
static std::string const input
Definition: EdmProvDump.cc:50
GlobalPoint getPosition(const DetId &id) const
Definition: RecHitTools.cc:140
T sqrt(T t)
Definition: SSEVec.h:23
bool isSilicon(const DetId &) const
Definition: RecHitTools.cc:478
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double f[11][100]
int findAndAssignTracksters(const TILES &, const std::vector< std::pair< int, int >> &)
void dumpClusters(const TILES &tiles, const std::vector< std::pair< int, int >> &layerIdx2layerandSoa, const int) const
GlobalPoint getPositionLayer(int layer, bool nose=false) const
Definition: RecHitTools.cc:152
void setGeometry(CaloGeometry const &)
Definition: RecHitTools.cc:79
void calculateLocalDensity(const TILES &, const int layerId, const std::vector< std::pair< int, int >> &)
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
unsigned int lastLayerEE(bool nose=false) const
Definition: RecHitTools.h:76
void dumpTracksters(const std::vector< std::pair< int, int >> &layerIdx2layerandSoa, const int, const std::vector< Trackster > &) const
unsigned int lastLayer(bool nose=false) const
Definition: RecHitTools.h:80
unsigned int getLayerWithOffset(const DetId &) const
Definition: RecHitTools.cc:381

◆ reset()

template<typename TILES >
void ticl::PatternRecognitionbyCLUE3D< TILES >::reset ( void  )
inlineprivate

Definition at line 96 of file PatternRecognitionbyCLUE3D.h.

References DummyCfis::c, and ticl::PatternRecognitionbyCLUE3D< TILES >::clusters_.

96  {
97  for (auto& c : clusters_) {
98  c.clear();
99  c.shrink_to_fit();
100  }
101  }
std::vector< ClustersOnLayer > clusters_

Member Data Documentation

◆ caloGeomToken_

template<typename TILES >
edm::ESGetToken<CaloGeometry, CaloGeometryRecord> ticl::PatternRecognitionbyCLUE3D< TILES >::caloGeomToken_
private

Definition at line 117 of file PatternRecognitionbyCLUE3D.h.

◆ clusters_

template<typename TILES >
std::vector<ClustersOnLayer> ticl::PatternRecognitionbyCLUE3D< TILES >::clusters_
private

◆ computeLocalTime_

template<typename TILES >
const bool ticl::PatternRecognitionbyCLUE3D< TILES >::computeLocalTime_
private

Definition at line 143 of file PatternRecognitionbyCLUE3D.h.

◆ criticalDensity_

template<typename TILES >
const std::vector<double> ticl::PatternRecognitionbyCLUE3D< TILES >::criticalDensity_
private

Definition at line 118 of file PatternRecognitionbyCLUE3D.h.

◆ criticalEtaPhiDistance_

template<typename TILES >
const std::vector<double> ticl::PatternRecognitionbyCLUE3D< TILES >::criticalEtaPhiDistance_
private

Definition at line 129 of file PatternRecognitionbyCLUE3D.h.

◆ criticalSelfDensity_

template<typename TILES >
const std::vector<double> ticl::PatternRecognitionbyCLUE3D< TILES >::criticalSelfDensity_
private

Definition at line 119 of file PatternRecognitionbyCLUE3D.h.

◆ criticalXYDistance_

template<typename TILES >
const std::vector<double> ticl::PatternRecognitionbyCLUE3D< TILES >::criticalXYDistance_
private

Definition at line 130 of file PatternRecognitionbyCLUE3D.h.

◆ criticalZDistanceLyr_

template<typename TILES >
const std::vector<int> ticl::PatternRecognitionbyCLUE3D< TILES >::criticalZDistanceLyr_
private

Definition at line 131 of file PatternRecognitionbyCLUE3D.h.

◆ cutHadProb_

template<typename TILES >
const float ticl::PatternRecognitionbyCLUE3D< TILES >::cutHadProb_
private

Definition at line 135 of file PatternRecognitionbyCLUE3D.h.

◆ densityEtaPhiDistanceSqr_

template<typename TILES >
const std::vector<double> ticl::PatternRecognitionbyCLUE3D< TILES >::densityEtaPhiDistanceSqr_
private

Definition at line 121 of file PatternRecognitionbyCLUE3D.h.

◆ densityOnSameLayer_

template<typename TILES >
const bool ticl::PatternRecognitionbyCLUE3D< TILES >::densityOnSameLayer_
private

Definition at line 124 of file PatternRecognitionbyCLUE3D.h.

◆ densitySiblingLayers_

template<typename TILES >
const std::vector<int> ticl::PatternRecognitionbyCLUE3D< TILES >::densitySiblingLayers_
private

Definition at line 120 of file PatternRecognitionbyCLUE3D.h.

◆ densityXYDistanceSqr_

template<typename TILES >
const std::vector<double> ticl::PatternRecognitionbyCLUE3D< TILES >::densityXYDistanceSqr_
private

Definition at line 122 of file PatternRecognitionbyCLUE3D.h.

◆ doPidCut_

template<typename TILES >
const bool ticl::PatternRecognitionbyCLUE3D< TILES >::doPidCut_
private

Definition at line 134 of file PatternRecognitionbyCLUE3D.h.

◆ eidInputName_

template<typename TILES >
const std::string ticl::PatternRecognitionbyCLUE3D< TILES >::eidInputName_
private

Definition at line 137 of file PatternRecognitionbyCLUE3D.h.

◆ eidMinClusterEnergy_

template<typename TILES >
const float ticl::PatternRecognitionbyCLUE3D< TILES >::eidMinClusterEnergy_
private

Definition at line 140 of file PatternRecognitionbyCLUE3D.h.

◆ eidNClusters_

template<typename TILES >
const int ticl::PatternRecognitionbyCLUE3D< TILES >::eidNClusters_
private

Definition at line 142 of file PatternRecognitionbyCLUE3D.h.

◆ eidNFeatures_

template<typename TILES >
const int ticl::PatternRecognitionbyCLUE3D< TILES >::eidNFeatures_ = 3
staticprivate

Definition at line 149 of file PatternRecognitionbyCLUE3D.h.

◆ eidNLayers_

template<typename TILES >
const int ticl::PatternRecognitionbyCLUE3D< TILES >::eidNLayers_
private

Definition at line 141 of file PatternRecognitionbyCLUE3D.h.

◆ eidOutputNameEnergy_

template<typename TILES >
const std::string ticl::PatternRecognitionbyCLUE3D< TILES >::eidOutputNameEnergy_
private

Definition at line 138 of file PatternRecognitionbyCLUE3D.h.

◆ eidOutputNameId_

template<typename TILES >
const std::string ticl::PatternRecognitionbyCLUE3D< TILES >::eidOutputNameId_
private

Definition at line 139 of file PatternRecognitionbyCLUE3D.h.

◆ eidSession_

template<typename TILES >
tensorflow::Session* ticl::PatternRecognitionbyCLUE3D< TILES >::eidSession_
private

Definition at line 147 of file PatternRecognitionbyCLUE3D.h.

◆ filter_on_categories_

template<typename TILES >
const std::vector<int> ticl::PatternRecognitionbyCLUE3D< TILES >::filter_on_categories_
private

Definition at line 136 of file PatternRecognitionbyCLUE3D.h.

◆ kernelDensityFactor_

template<typename TILES >
const std::vector<double> ticl::PatternRecognitionbyCLUE3D< TILES >::kernelDensityFactor_
private

Definition at line 123 of file PatternRecognitionbyCLUE3D.h.

◆ layersPosZ_

template<typename TILES >
std::vector<float> ticl::PatternRecognitionbyCLUE3D< TILES >::layersPosZ_
private

Definition at line 114 of file PatternRecognitionbyCLUE3D.h.

◆ minNumLayerCluster_

template<typename TILES >
const std::vector<int> ticl::PatternRecognitionbyCLUE3D< TILES >::minNumLayerCluster_
private

Definition at line 133 of file PatternRecognitionbyCLUE3D.h.

◆ nearestHigherOnSameLayer_

template<typename TILES >
const bool ticl::PatternRecognitionbyCLUE3D< TILES >::nearestHigherOnSameLayer_
private

Definition at line 125 of file PatternRecognitionbyCLUE3D.h.

◆ outlierMultiplier_

template<typename TILES >
const std::vector<double> ticl::PatternRecognitionbyCLUE3D< TILES >::outlierMultiplier_
private

Definition at line 132 of file PatternRecognitionbyCLUE3D.h.

◆ rescaleDensityByZ_

template<typename TILES >
const bool ticl::PatternRecognitionbyCLUE3D< TILES >::rescaleDensityByZ_
private

Definition at line 128 of file PatternRecognitionbyCLUE3D.h.

◆ rhtools_

template<typename TILES >
hgcal::RecHitTools ticl::PatternRecognitionbyCLUE3D< TILES >::rhtools_
private

Definition at line 146 of file PatternRecognitionbyCLUE3D.h.

◆ tracksterSeedAlgoId_

template<typename TILES >
std::vector<int> ticl::PatternRecognitionbyCLUE3D< TILES >::tracksterSeedAlgoId_
private

Definition at line 115 of file PatternRecognitionbyCLUE3D.h.

◆ useAbsoluteProjectiveScale_

template<typename TILES >
const bool ticl::PatternRecognitionbyCLUE3D< TILES >::useAbsoluteProjectiveScale_
private

Definition at line 126 of file PatternRecognitionbyCLUE3D.h.

◆ useClusterDimensionXY_

template<typename TILES >
const bool ticl::PatternRecognitionbyCLUE3D< TILES >::useClusterDimensionXY_
private

Definition at line 127 of file PatternRecognitionbyCLUE3D.h.

◆ usePCACleaning_

template<typename TILES >
const bool ticl::PatternRecognitionbyCLUE3D< TILES >::usePCACleaning_
private

Definition at line 144 of file PatternRecognitionbyCLUE3D.h.