File indexing completed on 2023-03-17 11:23:53
0001
0002
0003
0004 #include "LCToCPAssociatorByEnergyScoreImpl.h"
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0008 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0009
0010 #include "SimCalorimetry/HGCalAssociatorProducers/interface/AssociatorTools.h"
0011
0012 LCToCPAssociatorByEnergyScoreImpl::LCToCPAssociatorByEnergyScoreImpl(
0013 edm::EDProductGetter const& productGetter,
0014 bool hardScatterOnly,
0015 std::shared_ptr<hgcal::RecHitTools> recHitTools,
0016 const std::unordered_map<DetId, const HGCRecHit*>* hitMap)
0017 : hardScatterOnly_(hardScatterOnly), recHitTools_(recHitTools), hitMap_(hitMap), productGetter_(&productGetter) {
0018 layers_ = recHitTools_->lastLayerBH();
0019 }
0020
0021 hgcal::association LCToCPAssociatorByEnergyScoreImpl::makeConnections(
0022 const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0023
0024 const auto& clusters = *cCCH.product();
0025 const auto& caloParticles = *cPCH.product();
0026 auto nLayerClusters = clusters.size();
0027
0028
0029 std::vector<size_t> cPIndices;
0030 removeCPFromPU(caloParticles, cPIndices, hardScatterOnly_);
0031 auto nCaloParticles = cPIndices.size();
0032
0033
0034
0035
0036
0037 hgcal::caloParticleToLayerCluster cPOnLayer;
0038 cPOnLayer.resize(nCaloParticles);
0039 for (unsigned int i = 0; i < nCaloParticles; ++i) {
0040 cPOnLayer[i].resize(layers_ * 2);
0041 for (unsigned int j = 0; j < layers_ * 2; ++j) {
0042 cPOnLayer[i][j].caloParticleId = i;
0043 cPOnLayer[i][j].energy = 0.f;
0044 cPOnLayer[i][j].hits_and_fractions.clear();
0045
0046 }
0047 }
0048
0049
0050
0051
0052
0053
0054 std::unordered_map<DetId, std::vector<hgcal::detIdInfoInCluster>> detIdToCaloParticleId_Map;
0055 for (const auto& cpId : cPIndices) {
0056 const SimClusterRefVector& simClusterRefVector = caloParticles[cpId].simClusters();
0057 for (const auto& it_sc : simClusterRefVector) {
0058 const SimCluster& simCluster = (*(it_sc));
0059 const auto& hits_and_fractions = simCluster.hits_and_fractions();
0060 for (const auto& it_haf : hits_and_fractions) {
0061 const auto hitid = (it_haf.first);
0062 const auto cpLayerId =
0063 recHitTools_->getLayerWithOffset(hitid) + layers_ * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
0064 const auto itcheck = hitMap_->find(hitid);
0065 if (itcheck != hitMap_->end()) {
0066 auto hit_find_it = detIdToCaloParticleId_Map.find(hitid);
0067 if (hit_find_it == detIdToCaloParticleId_Map.end()) {
0068 detIdToCaloParticleId_Map[hitid] = std::vector<hgcal::detIdInfoInCluster>();
0069 detIdToCaloParticleId_Map[hitid].emplace_back(cpId, it_haf.second);
0070 } else {
0071 auto findHitIt = std::find(detIdToCaloParticleId_Map[hitid].begin(),
0072 detIdToCaloParticleId_Map[hitid].end(),
0073 hgcal::detIdInfoInCluster{cpId, it_haf.second});
0074 if (findHitIt != detIdToCaloParticleId_Map[hitid].end()) {
0075 findHitIt->fraction += it_haf.second;
0076 } else {
0077 detIdToCaloParticleId_Map[hitid].emplace_back(cpId, it_haf.second);
0078 }
0079 }
0080 const HGCRecHit* hit = itcheck->second;
0081 cPOnLayer[cpId][cpLayerId].energy += it_haf.second * hit->energy();
0082
0083
0084
0085
0086
0087
0088 auto& haf = cPOnLayer[cpId][cpLayerId].hits_and_fractions;
0089 auto found = std::find_if(
0090 std::begin(haf), std::end(haf), [&hitid](const std::pair<DetId, float>& v) { return v.first == hitid; });
0091 if (found != haf.end()) {
0092 found->second += it_haf.second;
0093 } else {
0094 cPOnLayer[cpId][cpLayerId].hits_and_fractions.emplace_back(hitid, it_haf.second);
0095 }
0096 }
0097 }
0098 }
0099 }
0100
0101 #ifdef EDM_ML_DEBUG
0102 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "cPOnLayer INFO" << std::endl;
0103 for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) {
0104 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "For CaloParticle Idx: " << cp << " we have: " << std::endl;
0105 for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) {
0106 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << " On Layer: " << cpp << " we have:" << std::endl;
0107 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0108 << " CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl;
0109 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0110 << " Energy: " << cPOnLayer[cp][cpp].energy << std::endl;
0111 double tot_energy = 0.;
0112 for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) {
0113 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0114 << " Hits/fraction/energy: " << (uint32_t)haf.first << "/" << haf.second << "/"
0115 << haf.second * hitMap_->at(haf.first)->energy() << std::endl;
0116 tot_energy += haf.second * hitMap_->at(haf.first)->energy();
0117 }
0118 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << " Tot Sum haf: " << tot_energy << std::endl;
0119 for (auto const& lc : cPOnLayer[cp][cpp].layerClusterIdToEnergyAndScore) {
0120 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << " lcIdx/energy/score: " << lc.first << "/"
0121 << lc.second.first << "/" << lc.second.second << std::endl;
0122 }
0123 }
0124 }
0125
0126 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "detIdToCaloParticleId_Map INFO" << std::endl;
0127 for (auto const& cp : detIdToCaloParticleId_Map) {
0128 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0129 << "For detId: " << (uint32_t)cp.first
0130 << " we have found the following connections with CaloParticles:" << std::endl;
0131 for (auto const& cpp : cp.second) {
0132 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0133 << " CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction
0134 << " and energy: " << cpp.fraction * hitMap_->at(cp.first)->energy() << std::endl;
0135 }
0136 }
0137 #endif
0138
0139
0140 std::unordered_map<DetId, std::vector<hgcal::detIdInfoInCluster>> detIdToLayerClusterId_Map;
0141
0142
0143
0144
0145 hgcal::layerClusterToCaloParticle cpsInLayerCluster;
0146 cpsInLayerCluster.resize(nLayerClusters);
0147
0148 for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0149 const std::vector<std::pair<DetId, float>>& hits_and_fractions = clusters[lcId].hitsAndFractions();
0150 unsigned int numberOfHitsInLC = hits_and_fractions.size();
0151 const auto firstHitDetId = hits_and_fractions[0].first;
0152 int lcLayerId =
0153 recHitTools_->getLayerWithOffset(firstHitDetId) + layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
0154
0155 for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0156 const auto rh_detid = hits_and_fractions[hitId].first;
0157 const auto rhFraction = hits_and_fractions[hitId].second;
0158
0159 auto hit_find_in_LC = detIdToLayerClusterId_Map.find(rh_detid);
0160 if (hit_find_in_LC == detIdToLayerClusterId_Map.end()) {
0161 detIdToLayerClusterId_Map[rh_detid] = std::vector<hgcal::detIdInfoInCluster>();
0162 }
0163 detIdToLayerClusterId_Map[rh_detid].emplace_back(lcId, rhFraction);
0164
0165 auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
0166
0167 if (hit_find_in_CP != detIdToCaloParticleId_Map.end()) {
0168 const auto itcheck = hitMap_->find(rh_detid);
0169 const HGCRecHit* hit = itcheck->second;
0170 for (auto& h : hit_find_in_CP->second) {
0171 cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[lcId].first += h.fraction * hit->energy();
0172 cpsInLayerCluster[lcId].emplace_back(h.clusterId, 0.f);
0173 }
0174 }
0175 }
0176 }
0177
0178 #ifdef EDM_ML_DEBUG
0179 for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0180 const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
0181 unsigned int numberOfHitsInLC = hits_and_fractions.size();
0182 const auto firstHitDetId = hits_and_fractions[0].first;
0183 int lcLayerId =
0184 recHitTools_->getLayerWithOffset(firstHitDetId) + layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194 std::vector<int> hitsToCaloParticleId(numberOfHitsInLC);
0195
0196 int maxCPId_byNumberOfHits = -1;
0197
0198 unsigned int maxCPNumberOfHitsInLC = 0;
0199
0200 int maxCPId_byEnergy = -1;
0201
0202 float maxEnergySharedLCandCP = 0.f;
0203
0204 float energyFractionOfLCinCP = 0.f;
0205
0206 float energyFractionOfCPinLC = 0.f;
0207 std::unordered_map<unsigned, unsigned> occurrencesCPinLC;
0208 unsigned int numberOfNoiseHitsInLC = 0;
0209 std::unordered_map<unsigned, float> CPEnergyInLC;
0210
0211 for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0212 const auto rh_detid = hits_and_fractions[hitId].first;
0213 const auto rhFraction = hits_and_fractions[hitId].second;
0214
0215 auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
0216
0217
0218
0219
0220
0221
0222
0223 if (rhFraction == 0.) {
0224 hitsToCaloParticleId[hitId] = -2;
0225 }
0226 if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) {
0227 hitsToCaloParticleId[hitId] -= 1;
0228 } else {
0229 const auto itcheck = hitMap_->find(rh_detid);
0230 const HGCRecHit* hit = itcheck->second;
0231 auto maxCPEnergyInLC = 0.f;
0232 auto maxCPId = -1;
0233 for (auto& h : hit_find_in_CP->second) {
0234 CPEnergyInLC[h.clusterId] += h.fraction * hit->energy();
0235
0236
0237 if (CPEnergyInLC[h.clusterId] > maxCPEnergyInLC) {
0238 maxCPEnergyInLC = CPEnergyInLC[h.clusterId];
0239 maxCPId = h.clusterId;
0240 }
0241 }
0242 hitsToCaloParticleId[hitId] = maxCPId;
0243 }
0244 }
0245
0246 for (const auto& c : hitsToCaloParticleId) {
0247 if (c < 0) {
0248 numberOfNoiseHitsInLC++;
0249 } else {
0250 occurrencesCPinLC[c]++;
0251 }
0252 }
0253
0254 for (const auto& c : occurrencesCPinLC) {
0255 if (c.second > maxCPNumberOfHitsInLC) {
0256 maxCPId_byNumberOfHits = c.first;
0257 maxCPNumberOfHitsInLC = c.second;
0258 }
0259 }
0260
0261 for (const auto& c : CPEnergyInLC) {
0262 if (c.second > maxEnergySharedLCandCP) {
0263 maxCPId_byEnergy = c.first;
0264 maxEnergySharedLCandCP = c.second;
0265 }
0266 }
0267
0268 float totalCPEnergyOnLayer = 0.f;
0269 if (maxCPId_byEnergy >= 0) {
0270 totalCPEnergyOnLayer = cPOnLayer[maxCPId_byEnergy][lcLayerId].energy;
0271 energyFractionOfCPinLC = maxEnergySharedLCandCP / totalCPEnergyOnLayer;
0272 if (clusters[lcId].energy() > 0.f) {
0273 energyFractionOfLCinCP = maxEnergySharedLCandCP / clusters[lcId].energy();
0274 }
0275 }
0276
0277 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0278 << std::setw(10) << "LayerId:\t" << std::setw(12) << "layerCluster\t" << std::setw(10) << "lc energy\t"
0279 << std::setw(5) << "nhits\t" << std::setw(12) << "noise hits\t" << std::setw(22) << "maxCPId_byNumberOfHits\t"
0280 << std::setw(8) << "nhitsCP\t" << std::setw(13) << "maxCPId_byEnergy\t" << std::setw(20)
0281 << "maxEnergySharedLCandCP\t" << std::setw(22) << "totalCPEnergyOnLayer\t" << std::setw(22)
0282 << "energyFractionOfLCinCP\t" << std::setw(25) << "energyFractionOfCPinLC\t"
0283 << "\n";
0284 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0285 << std::setw(10) << lcLayerId << "\t" << std::setw(12) << lcId << "\t" << std::setw(10)
0286 << clusters[lcId].energy() << "\t" << std::setw(5) << numberOfHitsInLC << "\t" << std::setw(12)
0287 << numberOfNoiseHitsInLC << "\t" << std::setw(22) << maxCPId_byNumberOfHits << "\t" << std::setw(8)
0288 << maxCPNumberOfHitsInLC << "\t" << std::setw(13) << maxCPId_byEnergy << "\t" << std::setw(20)
0289 << maxEnergySharedLCandCP << "\t" << std::setw(22) << totalCPEnergyOnLayer << "\t" << std::setw(22)
0290 << energyFractionOfLCinCP << "\t" << std::setw(25) << energyFractionOfCPinLC << "\n";
0291 }
0292
0293 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "Improved cPOnLayer INFO" << std::endl;
0294 for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) {
0295 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "For CaloParticle Idx: " << cp << " we have: " << std::endl;
0296 for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) {
0297 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << " On Layer: " << cpp << " we have:" << std::endl;
0298 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0299 << " CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl;
0300 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0301 << " Energy: " << cPOnLayer[cp][cpp].energy << std::endl;
0302 double tot_energy = 0.;
0303 for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) {
0304 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0305 << " Hits/fraction/energy: " << (uint32_t)haf.first << "/" << haf.second << "/"
0306 << haf.second * hitMap_->at(haf.first)->energy() << std::endl;
0307 tot_energy += haf.second * hitMap_->at(haf.first)->energy();
0308 }
0309 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << " Tot Sum haf: " << tot_energy << std::endl;
0310 for (auto const& lc : cPOnLayer[cp][cpp].layerClusterIdToEnergyAndScore) {
0311 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << " lcIdx/energy/score: " << lc.first << "/"
0312 << lc.second.first << "/" << lc.second.second << std::endl;
0313 }
0314 }
0315 }
0316
0317 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "Improved detIdToCaloParticleId_Map INFO" << std::endl;
0318 for (auto const& cp : detIdToCaloParticleId_Map) {
0319 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0320 << "For detId: " << (uint32_t)cp.first
0321 << " we have found the following connections with CaloParticles:" << std::endl;
0322 for (auto const& cpp : cp.second) {
0323 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0324 << " CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction
0325 << " and energy: " << cpp.fraction * hitMap_->at(cp.first)->energy() << std::endl;
0326 }
0327 }
0328 #endif
0329
0330
0331
0332 for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0333
0334 std::sort(cpsInLayerCluster[lcId].begin(), cpsInLayerCluster[lcId].end());
0335 auto last = std::unique(cpsInLayerCluster[lcId].begin(), cpsInLayerCluster[lcId].end());
0336 cpsInLayerCluster[lcId].erase(last, cpsInLayerCluster[lcId].end());
0337 const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
0338 unsigned int numberOfHitsInLC = hits_and_fractions.size();
0339
0340
0341 if (clusters[lcId].energy() == 0. && !cpsInLayerCluster[lcId].empty()) {
0342 for (auto& cpPair : cpsInLayerCluster[lcId]) {
0343 cpPair.second = 1.;
0344 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "layerClusterId : \t " << lcId << "\t CP id : \t"
0345 << cpPair.first << "\t score \t " << cpPair.second << "\n";
0346 }
0347 continue;
0348 }
0349
0350
0351
0352 float invLayerClusterEnergyWeight = 0.f;
0353 for (auto const& haf : hits_and_fractions) {
0354 invLayerClusterEnergyWeight +=
0355 (haf.second * hitMap_->at(haf.first)->energy()) * (haf.second * hitMap_->at(haf.first)->energy());
0356 }
0357 invLayerClusterEnergyWeight = 1.f / invLayerClusterEnergyWeight;
0358 for (unsigned int i = 0; i < numberOfHitsInLC; ++i) {
0359 DetId rh_detid = hits_and_fractions[i].first;
0360 float rhFraction = hits_and_fractions[i].second;
0361
0362 bool hitWithNoCP = (detIdToCaloParticleId_Map.find(rh_detid) == detIdToCaloParticleId_Map.end());
0363
0364 auto itcheck = hitMap_->find(rh_detid);
0365 const HGCRecHit* hit = itcheck->second;
0366 float hitEnergyWeight = hit->energy() * hit->energy();
0367
0368 for (auto& cpPair : cpsInLayerCluster[lcId]) {
0369 float cpFraction = 0.f;
0370 if (!hitWithNoCP) {
0371 auto findHitIt = std::find(detIdToCaloParticleId_Map[rh_detid].begin(),
0372 detIdToCaloParticleId_Map[rh_detid].end(),
0373 hgcal::detIdInfoInCluster{cpPair.first, 0.f});
0374 if (findHitIt != detIdToCaloParticleId_Map[rh_detid].end())
0375 cpFraction = findHitIt->fraction;
0376 }
0377 cpPair.second += std::min(std::pow(rhFraction - cpFraction, 2), std::pow(rhFraction, 2)) * hitEnergyWeight *
0378 invLayerClusterEnergyWeight;
0379 }
0380 }
0381 #ifdef EDM_ML_DEBUG
0382 if (cpsInLayerCluster[lcId].empty())
0383 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "layerCluster Id: \t" << lcId << "\tCP id:\t-1 "
0384 << "\t score \t-1\n";
0385 #endif
0386 }
0387
0388
0389 for (const auto& cpId : cPIndices) {
0390 for (unsigned int layerId = 0; layerId < layers_ * 2; ++layerId) {
0391 unsigned int CPNumberOfHits = cPOnLayer[cpId][layerId].hits_and_fractions.size();
0392 if (CPNumberOfHits == 0)
0393 continue;
0394 #ifdef EDM_ML_DEBUG
0395 int lcWithMaxEnergyInCP = -1;
0396 float maxEnergyLCinCP = 0.f;
0397 float CPenergy = cPOnLayer[cpId][layerId].energy;
0398 float CPEnergyFractionInLC = 0.f;
0399 for (auto& lc : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
0400 if (lc.second.first > maxEnergyLCinCP) {
0401 maxEnergyLCinCP = lc.second.first;
0402 lcWithMaxEnergyInCP = lc.first;
0403 }
0404 }
0405 if (CPenergy > 0.f)
0406 CPEnergyFractionInLC = maxEnergyLCinCP / CPenergy;
0407
0408 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0409 << std::setw(8) << "LayerId:\t" << std::setw(12) << "caloparticle\t" << std::setw(15) << "cp total energy\t"
0410 << std::setw(15) << "cpEnergyOnLayer\t" << std::setw(14) << "CPNhitsOnLayer\t" << std::setw(18)
0411 << "lcWithMaxEnergyInCP\t" << std::setw(15) << "maxEnergyLCinCP\t" << std::setw(20) << "CPEnergyFractionInLC"
0412 << "\n";
0413 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0414 << std::setw(8) << layerId << "\t" << std::setw(12) << cpId << "\t" << std::setw(15)
0415 << caloParticles[cpId].energy() << "\t" << std::setw(15) << CPenergy << "\t" << std::setw(14)
0416 << CPNumberOfHits << "\t" << std::setw(18) << lcWithMaxEnergyInCP << "\t" << std::setw(15) << maxEnergyLCinCP
0417 << "\t" << std::setw(20) << CPEnergyFractionInLC << "\n";
0418 #endif
0419
0420 float invCPEnergyWeight = 0.f;
0421 for (auto const& haf : cPOnLayer[cpId][layerId].hits_and_fractions) {
0422 invCPEnergyWeight += std::pow(haf.second * hitMap_->at(haf.first)->energy(), 2);
0423 }
0424 invCPEnergyWeight = 1.f / invCPEnergyWeight;
0425 for (unsigned int i = 0; i < CPNumberOfHits; ++i) {
0426 auto& cp_hitDetId = cPOnLayer[cpId][layerId].hits_and_fractions[i].first;
0427 auto& cpFraction = cPOnLayer[cpId][layerId].hits_and_fractions[i].second;
0428
0429 bool hitWithNoLC = false;
0430 if (cpFraction == 0.f)
0431 continue;
0432 auto hit_find_in_LC = detIdToLayerClusterId_Map.find(cp_hitDetId);
0433 if (hit_find_in_LC == detIdToLayerClusterId_Map.end())
0434 hitWithNoLC = true;
0435 auto itcheck = hitMap_->find(cp_hitDetId);
0436 const HGCRecHit* hit = itcheck->second;
0437 float hitEnergyWeight = hit->energy() * hit->energy();
0438 for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
0439 unsigned int layerClusterId = lcPair.first;
0440 float lcFraction = 0.f;
0441
0442 if (!hitWithNoLC) {
0443 auto findHitIt = std::find(detIdToLayerClusterId_Map[cp_hitDetId].begin(),
0444 detIdToLayerClusterId_Map[cp_hitDetId].end(),
0445 hgcal::detIdInfoInCluster{layerClusterId, 0.f});
0446 if (findHitIt != detIdToLayerClusterId_Map[cp_hitDetId].end())
0447 lcFraction = findHitIt->fraction;
0448 }
0449 lcPair.second.second += std::min(std::pow(lcFraction - cpFraction, 2), std::pow(cpFraction, 2)) *
0450 hitEnergyWeight * invCPEnergyWeight;
0451 #ifdef EDM_ML_DEBUG
0452 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0453 << "cpDetId:\t" << (uint32_t)cp_hitDetId << "\tlayerClusterId:\t" << layerClusterId << "\t"
0454 << "lcfraction,cpfraction:\t" << lcFraction << ", " << cpFraction << "\t"
0455 << "hitEnergyWeight:\t" << hitEnergyWeight << "\t"
0456 << "current score:\t" << lcPair.second.second << "\t"
0457 << "invCPEnergyWeight:\t" << invCPEnergyWeight << "\n";
0458 #endif
0459 }
0460 }
0461 #ifdef EDM_ML_DEBUG
0462 if (cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore.empty())
0463 LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "CP Id: \t" << cpId << "\tLC id:\t-1 "
0464 << "\t score \t-1\n";
0465
0466 for (const auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
0467 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0468 << "CP Id: \t" << cpId << "\t LC id: \t" << lcPair.first << "\t score \t" << lcPair.second.second
0469 << "\t shared energy:\t" << lcPair.second.first << "\t shared energy fraction:\t"
0470 << (lcPair.second.first / CPenergy) << "\n";
0471 }
0472 #endif
0473 }
0474 }
0475
0476 return {cpsInLayerCluster, cPOnLayer};
0477 }
0478
0479 hgcal::RecoToSimCollection LCToCPAssociatorByEnergyScoreImpl::associateRecoToSim(
0480 const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0481 hgcal::RecoToSimCollection returnValue(productGetter_);
0482 const auto& links = makeConnections(cCCH, cPCH);
0483
0484 const auto& cpsInLayerCluster = std::get<0>(links);
0485 for (size_t lcId = 0; lcId < cpsInLayerCluster.size(); ++lcId) {
0486 for (auto& cpPair : cpsInLayerCluster[lcId]) {
0487 LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0488 << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first << "\t score \t" << cpPair.second << "\n";
0489
0490 returnValue.insert(edm::Ref<reco::CaloClusterCollection>(cCCH, lcId),
0491 std::make_pair(edm::Ref<CaloParticleCollection>(cPCH, cpPair.first),
0492 cpPair.second)
0493 );
0494 }
0495 }
0496 return returnValue;
0497 }
0498
0499 hgcal::SimToRecoCollection LCToCPAssociatorByEnergyScoreImpl::associateSimToReco(
0500 const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0501 hgcal::SimToRecoCollection returnValue(productGetter_);
0502 const auto& links = makeConnections(cCCH, cPCH);
0503 const auto& cPOnLayer = std::get<1>(links);
0504 for (size_t cpId = 0; cpId < cPOnLayer.size(); ++cpId) {
0505 for (size_t layerId = 0; layerId < cPOnLayer[cpId].size(); ++layerId) {
0506 for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
0507 returnValue.insert(
0508 edm::Ref<CaloParticleCollection>(cPCH, cpId),
0509 std::make_pair(edm::Ref<reco::CaloClusterCollection>(cCCH, lcPair.first),
0510 std::make_pair(lcPair.second.first, lcPair.second.second))
0511 );
0512 }
0513 }
0514 }
0515 return returnValue;
0516 }