CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
LayerClusterAssociatorByEnergyScoreImpl Class Reference

#include <LayerClusterAssociatorByEnergyScoreImpl.h>

Inheritance diagram for LayerClusterAssociatorByEnergyScoreImpl:
hgcal::LayerClusterToCaloParticleAssociatorBaseImpl

Public Member Functions

hgcal::RecoToSimCollection associateRecoToSim (const edm::Handle< reco::CaloClusterCollection > &cCH, const edm::Handle< CaloParticleCollection > &cPCH) const override
 Associate a LayerCluster to CaloParticles. More...
 
hgcal::SimToRecoCollection associateSimToReco (const edm::Handle< reco::CaloClusterCollection > &cCH, const edm::Handle< CaloParticleCollection > &cPCH) const override
 Associate a CaloParticle to LayerClusters. More...
 
 LayerClusterAssociatorByEnergyScoreImpl (edm::EDProductGetter const &, bool, std::shared_ptr< hgcal::RecHitTools >, const std::unordered_map< DetId, const HGCRecHit * > *&)
 
- Public Member Functions inherited from hgcal::LayerClusterToCaloParticleAssociatorBaseImpl
 LayerClusterToCaloParticleAssociatorBaseImpl ()
 Constructor. More...
 
virtual ~LayerClusterToCaloParticleAssociatorBaseImpl ()
 Destructor. More...
 

Private Member Functions

hgcal::association makeConnections (const edm::Handle< reco::CaloClusterCollection > &cCH, const edm::Handle< CaloParticleCollection > &cPCH) const
 

Private Attributes

const bool hardScatterOnly_
 
const std::unordered_map< DetId, const HGCRecHit * > * hitMap_
 
unsigned layers_
 
const edm::EDProductGetterproductGetter_
 
std::shared_ptr< hgcal::RecHitToolsrecHitTools_
 

Detailed Description

Definition at line 47 of file LayerClusterAssociatorByEnergyScoreImpl.h.

Constructor & Destructor Documentation

◆ LayerClusterAssociatorByEnergyScoreImpl()

LayerClusterAssociatorByEnergyScoreImpl::LayerClusterAssociatorByEnergyScoreImpl ( edm::EDProductGetter const &  productGetter,
bool  hardScatterOnly,
std::shared_ptr< hgcal::RecHitTools recHitTools,
const std::unordered_map< DetId, const HGCRecHit * > *&  hitMap 
)
explicit

Definition at line 12 of file LayerClusterAssociatorByEnergyScoreImpl.cc.

17  : hardScatterOnly_(hardScatterOnly), recHitTools_(recHitTools), hitMap_(hitMap), productGetter_(&productGetter) {
18  layers_ = recHitTools_->lastLayerBH();
19 }

References layers_, and recHitTools_.

Member Function Documentation

◆ associateRecoToSim()

hgcal::RecoToSimCollection LayerClusterAssociatorByEnergyScoreImpl::associateRecoToSim ( const edm::Handle< reco::CaloClusterCollection > &  cCH,
const edm::Handle< CaloParticleCollection > &  cPCH 
) const
overridevirtual

Associate a LayerCluster to CaloParticles.

Reimplemented from hgcal::LayerClusterToCaloParticleAssociatorBaseImpl.

Definition at line 481 of file LayerClusterAssociatorByEnergyScoreImpl.cc.

482  {
484  const auto& links = makeConnections(cCCH, cPCH);
485 
486  const auto& cpsInLayerCluster = std::get<0>(links);
487  for (size_t lcId = 0; lcId < cpsInLayerCluster.size(); ++lcId) {
488  for (auto& cpPair : cpsInLayerCluster[lcId]) {
489  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
490  << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first << "\t score \t" << cpPair.second << "\n";
491  // Fill AssociationMap
492  returnValue.insert(edm::Ref<reco::CaloClusterCollection>(cCCH, lcId), // Ref to LC
493  std::make_pair(edm::Ref<CaloParticleCollection>(cPCH, cpPair.first),
494  cpPair.second) // Pair <Ref to CP, score>
495  );
496  }
497  }
498  return returnValue;
499 }

References edm::AssociationMap< Tag >::insert(), electronStore::links, LogDebug, makeConnections(), and productGetter_.

◆ associateSimToReco()

hgcal::SimToRecoCollection LayerClusterAssociatorByEnergyScoreImpl::associateSimToReco ( const edm::Handle< reco::CaloClusterCollection > &  cCH,
const edm::Handle< CaloParticleCollection > &  cPCH 
) const
overridevirtual

Associate a CaloParticle to LayerClusters.

Reimplemented from hgcal::LayerClusterToCaloParticleAssociatorBaseImpl.

Definition at line 501 of file LayerClusterAssociatorByEnergyScoreImpl.cc.

502  {
504  const auto& links = makeConnections(cCCH, cPCH);
505  const auto& cPOnLayer = std::get<1>(links);
506  for (size_t cpId = 0; cpId < cPOnLayer.size(); ++cpId) {
507  for (size_t layerId = 0; layerId < cPOnLayer[cpId].size(); ++layerId) {
508  for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
509  returnValue.insert(
510  edm::Ref<CaloParticleCollection>(cPCH, cpId), // Ref to CP
511  std::make_pair(edm::Ref<reco::CaloClusterCollection>(cCCH, lcPair.first), // Pair <Ref to LC,
512  std::make_pair(lcPair.second.first, lcPair.second.second)) // pair <energy, score> >
513  );
514  }
515  }
516  }
517  return returnValue;
518 }

References edm::AssociationMap< Tag >::insert(), electronStore::links, makeConnections(), and productGetter_.

◆ makeConnections()

hgcal::association LayerClusterAssociatorByEnergyScoreImpl::makeConnections ( const edm::Handle< reco::CaloClusterCollection > &  cCH,
const edm::Handle< CaloParticleCollection > &  cPCH 
) const
private

Definition at line 21 of file LayerClusterAssociatorByEnergyScoreImpl.cc.

22  {
23  // 1. Extract collections and filter CaloParticles, if required
24  const auto& clusters = *cCCH.product();
25  const auto& caloParticles = *cPCH.product();
26  auto nLayerClusters = clusters.size();
27  //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution.
28  std::vector<size_t> cPIndices;
29  //Consider CaloParticles coming from the hard scatterer
30  //excluding the PU contribution and save the indices.
31  removeCPFromPU(caloParticles, cPIndices);
32  auto nCaloParticles = cPIndices.size();
33 
34  // Initialize cPOnLayer. To be returned outside, since it contains the
35  // information to compute the CaloParticle-To-LayerCluster score.
37  cPOnLayer.resize(nCaloParticles);
38  for (unsigned int i = 0; i < nCaloParticles; ++i) {
39  cPOnLayer[i].resize(layers_ * 2);
40  for (unsigned int j = 0; j < layers_ * 2; ++j) {
41  cPOnLayer[i][j].caloParticleId = i;
42  cPOnLayer[i][j].energy = 0.f;
43  cPOnLayer[i][j].hits_and_fractions.clear();
44  //cPOnLayer[i][j].layerClusterIdToEnergyAndScore.reserve(nLayerClusters); // Not necessary but may improve performance
45  }
46  }
47 
48  // Fill detIdToCaloParticleId_Map and update cPOnLayer
49  std::unordered_map<DetId, std::vector<hgcal::detIdInfoInCluster>> detIdToCaloParticleId_Map;
50  for (const auto& cpId : cPIndices) {
51  const SimClusterRefVector& simClusterRefVector = caloParticles[cpId].simClusters();
52  for (const auto& it_sc : simClusterRefVector) {
53  const SimCluster& simCluster = (*(it_sc));
54  const auto& hits_and_fractions = simCluster.hits_and_fractions();
55  for (const auto& it_haf : hits_and_fractions) {
56  const auto hitid = (it_haf.first);
57  const auto cpLayerId =
58  recHitTools_->getLayerWithOffset(hitid) + layers_ * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
59  const auto itcheck = hitMap_->find(hitid);
60  if (itcheck != hitMap_->end()) {
61  auto hit_find_it = detIdToCaloParticleId_Map.find(hitid);
62  if (hit_find_it == detIdToCaloParticleId_Map.end()) {
63  detIdToCaloParticleId_Map[hitid] = std::vector<hgcal::detIdInfoInCluster>();
64  detIdToCaloParticleId_Map[hitid].emplace_back(cpId, it_haf.second);
65  } else {
66  auto findHitIt = std::find(detIdToCaloParticleId_Map[hitid].begin(),
67  detIdToCaloParticleId_Map[hitid].end(),
68  hgcal::detIdInfoInCluster{cpId, it_haf.second});
69  if (findHitIt != detIdToCaloParticleId_Map[hitid].end()) {
70  findHitIt->fraction += it_haf.second;
71  } else {
72  detIdToCaloParticleId_Map[hitid].emplace_back(cpId, it_haf.second);
73  }
74  }
75  const HGCRecHit* hit = itcheck->second;
76  cPOnLayer[cpId][cpLayerId].energy += it_haf.second * hit->energy();
77  // We need to compress the hits and fractions in order to have a
78  // reasonable score between CP and LC. Imagine, for example, that a
79  // CP has detID X used by 2 SimClusters with different fractions. If
80  // a single LC uses X with fraction 1 and is compared to the 2
81  // contributions separately, it will be assigned a score != 0, which
82  // is wrong.
83  auto& haf = cPOnLayer[cpId][cpLayerId].hits_and_fractions;
84  auto found = std::find_if(
85  std::begin(haf), std::end(haf), [&hitid](const std::pair<DetId, float>& v) { return v.first == hitid; });
86  if (found != haf.end()) {
87  found->second += it_haf.second;
88  } else {
89  cPOnLayer[cpId][cpLayerId].hits_and_fractions.emplace_back(hitid, it_haf.second);
90  }
91  }
92  }
93  }
94  }
95 
96 #ifdef EDM_ML_DEBUG
97  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "cPOnLayer INFO" << std::endl;
98  for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) {
99  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "For CaloParticle Idx: " << cp << " we have: " << std::endl;
100  for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) {
101  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " On Layer: " << cpp << " we have:" << std::endl;
102  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
103  << " CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl;
104  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
105  << " Energy: " << cPOnLayer[cp][cpp].energy << std::endl;
106  double tot_energy = 0.;
107  for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) {
108  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
109  << " Hits/fraction/energy: " << (uint32_t)haf.first << "/" << haf.second << "/"
110  << haf.second * hitMap_->at(haf.first)->energy() << std::endl;
111  tot_energy += haf.second * hitMap_->at(haf.first)->energy();
112  }
113  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " Tot Sum haf: " << tot_energy << std::endl;
114  for (auto const& lc : cPOnLayer[cp][cpp].layerClusterIdToEnergyAndScore) {
115  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " lcIdx/energy/score: " << lc.first << "/"
116  << lc.second.first << "/" << lc.second.second << std::endl;
117  }
118  }
119  }
120 
121  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "detIdToCaloParticleId_Map INFO" << std::endl;
122  for (auto const& cp : detIdToCaloParticleId_Map) {
123  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
124  << "For detId: " << (uint32_t)cp.first
125  << " we have found the following connections with CaloParticles:" << std::endl;
126  for (auto const& cpp : cp.second) {
127  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
128  << " CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction
129  << " and energy: " << cpp.fraction * hitMap_->at(cp.first)->energy() << std::endl;
130  }
131  }
132 #endif
133 
134  // Fill detIdToLayerClusterId_Map and cpsInLayerCluster; update cPOnLayer
135  std::unordered_map<DetId, std::vector<hgcal::detIdInfoInCluster>> detIdToLayerClusterId_Map;
136  // this contains the ids of the caloparticles contributing with at least one
137  // hit to the layer cluster and the reconstruction error. To be returned
138  // since this contains the information to compute the
139  // LayerCluster-To-CaloParticle score.
140  hgcal::layerClusterToCaloParticle cpsInLayerCluster;
141  cpsInLayerCluster.resize(nLayerClusters);
142 
143  for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
144  const std::vector<std::pair<DetId, float>>& hits_and_fractions = clusters[lcId].hitsAndFractions();
145  unsigned int numberOfHitsInLC = hits_and_fractions.size();
146  const auto firstHitDetId = hits_and_fractions[0].first;
147  int lcLayerId =
148  recHitTools_->getLayerWithOffset(firstHitDetId) + layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
149 
150  for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
151  const auto rh_detid = hits_and_fractions[hitId].first;
152  const auto rhFraction = hits_and_fractions[hitId].second;
153 
154  auto hit_find_in_LC = detIdToLayerClusterId_Map.find(rh_detid);
155  if (hit_find_in_LC == detIdToLayerClusterId_Map.end()) {
156  detIdToLayerClusterId_Map[rh_detid] = std::vector<hgcal::detIdInfoInCluster>();
157  }
158  detIdToLayerClusterId_Map[rh_detid].emplace_back(lcId, rhFraction);
159 
160  auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
161 
162  if (hit_find_in_CP != detIdToCaloParticleId_Map.end()) {
163  const auto itcheck = hitMap_->find(rh_detid);
164  const HGCRecHit* hit = itcheck->second;
165  for (auto& h : hit_find_in_CP->second) {
166  cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[lcId].first += h.fraction * hit->energy();
167  cpsInLayerCluster[lcId].emplace_back(h.clusterId, 0.f);
168  }
169  }
170  } // End loop over hits on a LayerCluster
171  } // End of loop over LayerClusters
172 
173 #ifdef EDM_ML_DEBUG
174  for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
175  const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
176  unsigned int numberOfHitsInLC = hits_and_fractions.size();
177  const auto firstHitDetId = hits_and_fractions[0].first;
178  int lcLayerId =
179  recHitTools_->getLayerWithOffset(firstHitDetId) + layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
180 
181  // This vector will store, for each hit in the Layercluster, the index of
182  // the CaloParticle that contributed the most, in terms of energy, to it.
183  // Special values are:
184  //
185  // -2 --> the reconstruction fraction of the RecHit is 0 (used in the past to monitor Halo Hits)
186  // -3 --> same as before with the added condition that no CaloParticle has been linked to this RecHit
187  // -1 --> the reco fraction is >0, but no CaloParticle has been linked to it
188  // >=0 --> index of the linked CaloParticle
189  std::vector<int> hitsToCaloParticleId(numberOfHitsInLC);
190  // This will store the index of the CaloParticle linked to the LayerCluster that has the most number of hits in common.
191  int maxCPId_byNumberOfHits = -1;
192  // This will store the maximum number of shared hits between a Layercluster andd a CaloParticle
193  unsigned int maxCPNumberOfHitsInLC = 0;
194  // This will store the index of the CaloParticle linked to the LayerCluster that has the most energy in common.
195  int maxCPId_byEnergy = -1;
196  // This will store the maximum number of shared energy between a Layercluster and a CaloParticle
197  float maxEnergySharedLCandCP = 0.f;
198  // This will store the fraction of the LayerCluster energy shared with the best(energy) CaloParticle: e_shared/lc_energy
199  float energyFractionOfLCinCP = 0.f;
200  // This will store the fraction of the CaloParticle energy shared with the LayerCluster: e_shared/cp_energy
201  float energyFractionOfCPinLC = 0.f;
202  std::unordered_map<unsigned, unsigned> occurrencesCPinLC;
203  unsigned int numberOfNoiseHitsInLC = 0;
204  std::unordered_map<unsigned, float> CPEnergyInLC;
205 
206  for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
207  const auto rh_detid = hits_and_fractions[hitId].first;
208  const auto rhFraction = hits_and_fractions[hitId].second;
209 
210  auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
211 
212  // if the fraction is zero or the hit does not belong to any calo
213  // particle, set the caloparticleId for the hit to -1 this will
214  // contribute to the number of noise hits
215 
216  // MR Remove the case in which the fraction is 0, since this could be a
217  // real hit that has been marked as halo.
218  if (rhFraction == 0.) {
219  hitsToCaloParticleId[hitId] = -2;
220  }
221  if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) {
222  hitsToCaloParticleId[hitId] -= 1;
223  } else {
224  const auto itcheck = hitMap_->find(rh_detid);
225  const HGCRecHit* hit = itcheck->second;
226  auto maxCPEnergyInLC = 0.f;
227  auto maxCPId = -1;
228  for (auto& h : hit_find_in_CP->second) {
229  CPEnergyInLC[h.clusterId] += h.fraction * hit->energy();
230  // Keep track of which CaloParticle ccontributed the most, in terms
231  // of energy, to this specific LayerCluster.
232  if (CPEnergyInLC[h.clusterId] > maxCPEnergyInLC) {
233  maxCPEnergyInLC = CPEnergyInLC[h.clusterId];
234  maxCPId = h.clusterId;
235  }
236  }
237  hitsToCaloParticleId[hitId] = maxCPId;
238  }
239  } // End loop over hits on a LayerCluster
240 
241  for (const auto& c : hitsToCaloParticleId) {
242  if (c < 0) {
243  numberOfNoiseHitsInLC++;
244  } else {
245  occurrencesCPinLC[c]++;
246  }
247  }
248 
249  for (const auto& c : occurrencesCPinLC) {
250  if (c.second > maxCPNumberOfHitsInLC) {
251  maxCPId_byNumberOfHits = c.first;
252  maxCPNumberOfHitsInLC = c.second;
253  }
254  }
255 
256  for (const auto& c : CPEnergyInLC) {
257  if (c.second > maxEnergySharedLCandCP) {
258  maxCPId_byEnergy = c.first;
259  maxEnergySharedLCandCP = c.second;
260  }
261  }
262 
263  float totalCPEnergyOnLayer = 0.f;
264  if (maxCPId_byEnergy >= 0) {
265  totalCPEnergyOnLayer = cPOnLayer[maxCPId_byEnergy][lcLayerId].energy;
266  energyFractionOfCPinLC = maxEnergySharedLCandCP / totalCPEnergyOnLayer;
267  if (clusters[lcId].energy() > 0.f) {
268  energyFractionOfLCinCP = maxEnergySharedLCandCP / clusters[lcId].energy();
269  }
270  }
271 
272  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << std::setw(10) << "LayerId:"
273  << "\t" << std::setw(12) << "layerCluster"
274  << "\t" << std::setw(10) << "lc energy"
275  << "\t" << std::setw(5) << "nhits"
276  << "\t" << std::setw(12) << "noise hits"
277  << "\t" << std::setw(22) << "maxCPId_byNumberOfHits"
278  << "\t" << std::setw(8) << "nhitsCP"
279  << "\t" << std::setw(13) << "maxCPId_byEnergy"
280  << "\t" << std::setw(20) << "maxEnergySharedLCandCP"
281  << "\t" << std::setw(22) << "totalCPEnergyOnLayer"
282  << "\t" << std::setw(22) << "energyFractionOfLCinCP"
283  << "\t" << std::setw(25) << "energyFractionOfCPinLC"
284  << "\t"
285  << "\n";
286  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
287  << std::setw(10) << lcLayerId << "\t" << std::setw(12) << lcId << "\t" << std::setw(10)
288  << clusters[lcId].energy() << "\t" << std::setw(5) << numberOfHitsInLC << "\t" << std::setw(12)
289  << numberOfNoiseHitsInLC << "\t" << std::setw(22) << maxCPId_byNumberOfHits << "\t" << std::setw(8)
290  << maxCPNumberOfHitsInLC << "\t" << std::setw(13) << maxCPId_byEnergy << "\t" << std::setw(20)
291  << maxEnergySharedLCandCP << "\t" << std::setw(22) << totalCPEnergyOnLayer << "\t" << std::setw(22)
292  << energyFractionOfLCinCP << "\t" << std::setw(25) << energyFractionOfCPinLC << "\n";
293  } // End of loop over LayerClusters
294 
295  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "Improved cPOnLayer INFO" << std::endl;
296  for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) {
297  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "For CaloParticle Idx: " << cp << " we have: " << std::endl;
298  for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) {
299  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " On Layer: " << cpp << " we have:" << std::endl;
300  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
301  << " CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl;
302  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
303  << " Energy: " << cPOnLayer[cp][cpp].energy << std::endl;
304  double tot_energy = 0.;
305  for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) {
306  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
307  << " Hits/fraction/energy: " << (uint32_t)haf.first << "/" << haf.second << "/"
308  << haf.second * hitMap_->at(haf.first)->energy() << std::endl;
309  tot_energy += haf.second * hitMap_->at(haf.first)->energy();
310  }
311  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " Tot Sum haf: " << tot_energy << std::endl;
312  for (auto const& lc : cPOnLayer[cp][cpp].layerClusterIdToEnergyAndScore) {
313  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " lcIdx/energy/score: " << lc.first << "/"
314  << lc.second.first << "/" << lc.second.second << std::endl;
315  }
316  }
317  }
318 
319  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "Improved detIdToCaloParticleId_Map INFO" << std::endl;
320  for (auto const& cp : detIdToCaloParticleId_Map) {
321  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
322  << "For detId: " << (uint32_t)cp.first
323  << " we have found the following connections with CaloParticles:" << std::endl;
324  for (auto const& cpp : cp.second) {
325  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
326  << " CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction
327  << " and energy: " << cpp.fraction * hitMap_->at(cp.first)->energy() << std::endl;
328  }
329  }
330 #endif
331 
332  // Update cpsInLayerCluster; compute the score LayerCluster-to-CaloParticle,
333  // together with the returned AssociationMap
334  for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
335  // find the unique caloparticles id contributing to the layer clusters
336  std::sort(cpsInLayerCluster[lcId].begin(), cpsInLayerCluster[lcId].end());
337  auto last = std::unique(cpsInLayerCluster[lcId].begin(), cpsInLayerCluster[lcId].end());
338  cpsInLayerCluster[lcId].erase(last, cpsInLayerCluster[lcId].end());
339  const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
340  unsigned int numberOfHitsInLC = hits_and_fractions.size();
341  // If a reconstructed LayerCluster has energy 0 but is linked to a
342  // CaloParticle, assigned score 1
343  if (clusters[lcId].energy() == 0. && !cpsInLayerCluster[lcId].empty()) {
344  for (auto& cpPair : cpsInLayerCluster[lcId]) {
345  cpPair.second = 1.;
346  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "layerClusterId : \t " << lcId << "\t CP id : \t"
347  << cpPair.first << "\t score \t " << cpPair.second << "\n";
348  }
349  continue;
350  }
351 
352  // Compute the correct normalization
353  float invLayerClusterEnergyWeight = 0.f;
354  for (auto const& haf : clusters[lcId].hitsAndFractions()) {
355  invLayerClusterEnergyWeight +=
356  (haf.second * hitMap_->at(haf.first)->energy()) * (haf.second * hitMap_->at(haf.first)->energy());
357  }
358  invLayerClusterEnergyWeight = 1.f / invLayerClusterEnergyWeight;
359  for (unsigned int i = 0; i < numberOfHitsInLC; ++i) {
360  DetId rh_detid = hits_and_fractions[i].first;
361  float rhFraction = hits_and_fractions[i].second;
362 
363  bool hitWithNoCP = (detIdToCaloParticleId_Map.find(rh_detid) == detIdToCaloParticleId_Map.end());
364 
365  auto itcheck = hitMap_->find(rh_detid);
366  const HGCRecHit* hit = itcheck->second;
367  float hitEnergyWeight = hit->energy() * hit->energy();
368 
369  for (auto& cpPair : cpsInLayerCluster[lcId]) {
370  float cpFraction = 0.f;
371  if (!hitWithNoCP) {
372  auto findHitIt = std::find(detIdToCaloParticleId_Map[rh_detid].begin(),
373  detIdToCaloParticleId_Map[rh_detid].end(),
374  hgcal::detIdInfoInCluster{cpPair.first, 0.f});
375  if (findHitIt != detIdToCaloParticleId_Map[rh_detid].end())
376  cpFraction = findHitIt->fraction;
377  }
378  cpPair.second +=
379  (rhFraction - cpFraction) * (rhFraction - cpFraction) * hitEnergyWeight * invLayerClusterEnergyWeight;
380  }
381  } // End of loop over Hits within a LayerCluster
382 #ifdef EDM_ML_DEBUG
383  if (cpsInLayerCluster[lcId].empty())
384  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "layerCluster Id: \t" << lcId << "\tCP id:\t-1 "
385  << "\t score \t-1"
386  << "\n";
387 #endif
388  } // End of loop over LayerClusters
389 
390  // Compute the CaloParticle-To-LayerCluster score
391  for (const auto& cpId : cPIndices) {
392  for (unsigned int layerId = 0; layerId < layers_ * 2; ++layerId) {
393  unsigned int CPNumberOfHits = cPOnLayer[cpId][layerId].hits_and_fractions.size();
394  if (CPNumberOfHits == 0)
395  continue;
396 #ifdef EDM_ML_DEBUG
397  int lcWithMaxEnergyInCP = -1;
398  float maxEnergyLCinCP = 0.f;
399  float CPenergy = cPOnLayer[cpId][layerId].energy;
400  float CPEnergyFractionInLC = 0.f;
401  for (auto& lc : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
402  if (lc.second.first > maxEnergyLCinCP) {
403  maxEnergyLCinCP = lc.second.first;
404  lcWithMaxEnergyInCP = lc.first;
405  }
406  }
407  if (CPenergy > 0.f)
408  CPEnergyFractionInLC = maxEnergyLCinCP / CPenergy;
409 
410  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
411  << std::setw(8) << "LayerId:\t" << std::setw(12) << "caloparticle\t" << std::setw(15) << "cp total energy\t"
412  << std::setw(15) << "cpEnergyOnLayer\t" << std::setw(14) << "CPNhitsOnLayer\t" << std::setw(18)
413  << "lcWithMaxEnergyInCP\t" << std::setw(15) << "maxEnergyLCinCP\t" << std::setw(20) << "CPEnergyFractionInLC"
414  << "\n";
415  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
416  << std::setw(8) << layerId << "\t" << std::setw(12) << cpId << "\t" << std::setw(15)
417  << caloParticles[cpId].energy() << "\t" << std::setw(15) << CPenergy << "\t" << std::setw(14)
418  << CPNumberOfHits << "\t" << std::setw(18) << lcWithMaxEnergyInCP << "\t" << std::setw(15) << maxEnergyLCinCP
419  << "\t" << std::setw(20) << CPEnergyFractionInLC << "\n";
420 #endif
421  // Compute the correct normalization
422  float invCPEnergyWeight = 0.f;
423  for (auto const& haf : cPOnLayer[cpId][layerId].hits_and_fractions) {
424  invCPEnergyWeight += std::pow(haf.second * hitMap_->at(haf.first)->energy(), 2);
425  }
426  invCPEnergyWeight = 1.f / invCPEnergyWeight;
427  for (unsigned int i = 0; i < CPNumberOfHits; ++i) {
428  auto& cp_hitDetId = cPOnLayer[cpId][layerId].hits_and_fractions[i].first;
429  auto& cpFraction = cPOnLayer[cpId][layerId].hits_and_fractions[i].second;
430 
431  bool hitWithNoLC = false;
432  if (cpFraction == 0.f)
433  continue; //hopefully this should never happen
434  auto hit_find_in_LC = detIdToLayerClusterId_Map.find(cp_hitDetId);
435  if (hit_find_in_LC == detIdToLayerClusterId_Map.end())
436  hitWithNoLC = true;
437  auto itcheck = hitMap_->find(cp_hitDetId);
438  const HGCRecHit* hit = itcheck->second;
439  float hitEnergyWeight = hit->energy() * hit->energy();
440  for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
441  unsigned int layerClusterId = lcPair.first;
442  float lcFraction = 0.f;
443 
444  if (!hitWithNoLC) {
445  auto findHitIt = std::find(detIdToLayerClusterId_Map[cp_hitDetId].begin(),
446  detIdToLayerClusterId_Map[cp_hitDetId].end(),
447  hgcal::detIdInfoInCluster{layerClusterId, 0.f});
448  if (findHitIt != detIdToLayerClusterId_Map[cp_hitDetId].end())
449  lcFraction = findHitIt->fraction;
450  }
451  lcPair.second.second +=
452  (lcFraction - cpFraction) * (lcFraction - cpFraction) * hitEnergyWeight * invCPEnergyWeight;
453 #ifdef EDM_ML_DEBUG
454  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
455  << "cpDetId:\t" << (uint32_t)cp_hitDetId << "\tlayerClusterId:\t" << layerClusterId << "\t"
456  << "lcfraction,cpfraction:\t" << lcFraction << ", " << cpFraction << "\t"
457  << "hitEnergyWeight:\t" << hitEnergyWeight << "\t"
458  << "current score:\t" << lcPair.second.second << "\t"
459  << "invCPEnergyWeight:\t" << invCPEnergyWeight << "\n";
460 #endif
461  } // End of loop over LayerClusters linked to hits of this CaloParticle
462  } // End of loop over hits of CaloParticle on a Layer
463 #ifdef EDM_ML_DEBUG
464  if (cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore.empty())
465  LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "CP Id: \t" << cpId << "\tLC id:\t-1 "
466  << "\t score \t-1"
467  << "\n";
468 
469  for (const auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
470  LogDebug("LayerClusterAssociatorByEnergyScoreImpl")
471  << "CP Id: \t" << cpId << "\t LC id: \t" << lcPair.first << "\t score \t" << lcPair.second.second
472  << "\t shared energy:\t" << lcPair.second.first << "\t shared energy fraction:\t"
473  << (lcPair.second.first / CPenergy) << "\n";
474  }
475 #endif
476  }
477  }
478  return {cpsInLayerCluster, cPOnLayer};
479 }

References c, caloTruthCellsProducer_cfi::caloParticles, bsc_activity_cfg::clusters, CommonMethods::cp(), relativeConstraints::empty, mps_fire::end, HCALHighEnergyHPDFilter_cfi::energy, f, spr::find(), cms::cuda::for(), newFWLiteAna::found, hitMap_, SimCluster::hits_and_fractions(), mps_fire::i, dqmiolumiharvest::j, dqmdumpme::last, layers_, LogDebug, funct::pow(), edm::Handle< T >::product(), recHitTools_, removeCPFromPU(), jetUpdater_cfi::sort, tier0::unique(), and findQualityFiles::v.

Referenced by associateRecoToSim(), and associateSimToReco().

Member Data Documentation

◆ hardScatterOnly_

const bool LayerClusterAssociatorByEnergyScoreImpl::hardScatterOnly_
private

Definition at line 61 of file LayerClusterAssociatorByEnergyScoreImpl.h.

◆ hitMap_

const std::unordered_map<DetId, const HGCRecHit *>* LayerClusterAssociatorByEnergyScoreImpl::hitMap_
private

Definition at line 63 of file LayerClusterAssociatorByEnergyScoreImpl.h.

Referenced by makeConnections().

◆ layers_

unsigned LayerClusterAssociatorByEnergyScoreImpl::layers_
private

◆ productGetter_

const edm::EDProductGetter* LayerClusterAssociatorByEnergyScoreImpl::productGetter_
private

◆ recHitTools_

std::shared_ptr<hgcal::RecHitTools> LayerClusterAssociatorByEnergyScoreImpl::recHitTools_
private
removeCPFromPU
void removeCPFromPU(const std::vector< CaloParticle > &caloParticles, std::vector< size_t > &cPIndices)
Definition: AssociatorTools.h:7
mps_fire.i
i
Definition: mps_fire.py:428
edm::Handle::product
T const * product() const
Definition: Handle.h:70
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
LayerClusterAssociatorByEnergyScoreImpl::productGetter_
const edm::EDProductGetter * productGetter_
Definition: LayerClusterAssociatorByEnergyScoreImpl.h:65
hgcal_conditions::parameters
Definition: HGCConditions.h:86
edm::RefVector< SimClusterCollection >
SimCluster
Monte Carlo truth information used for tracking validation.
Definition: SimCluster.h:29
findQualityFiles.v
v
Definition: findQualityFiles.py:179
LayerClusterAssociatorByEnergyScoreImpl::hardScatterOnly_
const bool hardScatterOnly_
Definition: LayerClusterAssociatorByEnergyScoreImpl.h:61
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
newFWLiteAna.found
found
Definition: newFWLiteAna.py:118
LayerClusterAssociatorByEnergyScoreImpl::recHitTools_
std::shared_ptr< hgcal::RecHitTools > recHitTools_
Definition: LayerClusterAssociatorByEnergyScoreImpl.h:62
edm::Ref
Definition: AssociativeIterator.h:58
DetId
Definition: DetId.h:17
dqmdumpme.last
last
Definition: dqmdumpme.py:56
caloTruthCellsProducer_cfi.caloParticles
caloParticles
Definition: caloTruthCellsProducer_cfi.py:6
mps_fire.end
end
Definition: mps_fire.py:242
HCALHighEnergyHPDFilter_cfi.energy
energy
Definition: HCALHighEnergyHPDFilter_cfi.py:5
h
HGCRecHit
Definition: HGCRecHit.h:14
LayerClusterAssociatorByEnergyScoreImpl::layers_
unsigned layers_
Definition: LayerClusterAssociatorByEnergyScoreImpl.h:64
bsc_activity_cfg.clusters
clusters
Definition: bsc_activity_cfg.py:36
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
hgcal::detIdInfoInCluster
Definition: LayerClusterAssociatorByEnergyScoreImpl.h:18
edm::AssociationMap
Definition: AssociationMap.h:48
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
edm::refcoreimpl::productGetter
EDProductGetter const * productGetter(std::atomic< void const * > const &iCache)
Definition: refcore_implementation.h:72
cms::cuda::for
for(int i=first, nt=offsets[nh];i< nt;i+=gridDim.x *blockDim.x)
Definition: HistoContainer.h:27
tier0.unique
def unique(seq, keepstr=True)
Definition: tier0.py:24
LayerClusterAssociatorByEnergyScoreImpl::makeConnections
hgcal::association makeConnections(const edm::Handle< reco::CaloClusterCollection > &cCH, const edm::Handle< CaloParticleCollection > &cPCH) const
Definition: LayerClusterAssociatorByEnergyScoreImpl.cc:21
relativeConstraints.empty
bool empty
Definition: relativeConstraints.py:46
electronStore.links
links
Definition: electronStore.py:149
hgcal::caloParticleToLayerCluster
std::vector< std::vector< hgcal::caloParticleOnLayer > > caloParticleToLayerCluster
Definition: LayerClusterAssociatorByEnergyScoreImpl.h:43
SimCluster::hits_and_fractions
std::vector< std::pair< uint32_t, float > > hits_and_fractions() const
Returns list of rechit IDs and fractions for this SimCluster.
Definition: SimCluster.h:184
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
c
auto & c
Definition: CAHitNtupletGeneratorKernelsImpl.h:46
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
CommonMethods.cp
def cp(fromDir, toDir, listOfFiles, overwrite=False, smallList=False)
Definition: CommonMethods.py:192
LayerClusterAssociatorByEnergyScoreImpl::hitMap_
const std::unordered_map< DetId, const HGCRecHit * > * hitMap_
Definition: LayerClusterAssociatorByEnergyScoreImpl.h:63
hit
Definition: SiStripHitEffFromCalibTree.cc:88
hgcal::layerClusterToCaloParticle
std::vector< std::vector< std::pair< unsigned int, float > > > layerClusterToCaloParticle
Definition: LayerClusterAssociatorByEnergyScoreImpl.h:42