HGCalAssociatorProducers/plugins/LCToSCAssociatorByEnergyScoreImpl.cc

0001 #include "LCToSCAssociatorByEnergyScoreImpl.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0005
0006 template <typename HIT>
0007 LCToSCAssociatorByEnergyScoreImpl<HIT>::LCToSCAssociatorByEnergyScoreImpl(
0008     edm::EDProductGetter const& productGetter,
0009     bool hardScatterOnly,
0010     std::shared_ptr<hgcal::RecHitTools> recHitTools,
0011     const std::unordered_map<DetId, const unsigned int>* hitMap,
0012     std::vector<const HIT*>& hits)
0013     : hardScatterOnly_(hardScatterOnly),
0014       recHitTools_(recHitTools),
0015       hitMap_(hitMap),
0016       productGetter_(&productGetter),
0017       hits_(hits) {
0018   if constexpr (std::is_same_v<HIT, HGCRecHit>)
0019     layers_ = recHitTools_->lastLayerBH();
0020   else
0021     layers_ = 6;  //EB + 4 HB + HO
0022 }
0023
0024 template <typename HIT>
0025 ticl::association LCToSCAssociatorByEnergyScoreImpl<HIT>::makeConnections(
0026     const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<SimClusterCollection>& sCCH) const {
0027   // Get collections
0028   const auto& clusters = *cCCH.product();
0029   const auto& simClusters = *sCCH.product();
0030   auto nLayerClusters = clusters.size();
0031
0032   //There shouldn't be any SimTracks from different crossings, but maybe they will be added later.
0033   //At the moment there should be one SimTrack in each SimCluster.
0034   auto nSimClusters = simClusters.size();
0035   std::vector<size_t> sCIndices;
0036   for (unsigned int scId = 0; scId < nSimClusters; ++scId) {
0037     if (hardScatterOnly_ && (simClusters[scId].g4Tracks()[0].eventId().event() != 0 or
0038                              simClusters[scId].g4Tracks()[0].eventId().bunchCrossing() != 0)) {
0039       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0040           << "Excluding SimCluster from event: " << simClusters[scId].g4Tracks()[0].eventId().event()
0041           << " with BX: " << simClusters[scId].g4Tracks()[0].eventId().bunchCrossing() << std::endl;
0042       continue;
0043     }
0044     sCIndices.emplace_back(scId);
0045   }
0046   nSimClusters = sCIndices.size();
0047
0048   // Initialize lcsInSimCluster. It contains the simClusterOnLayer structure for all simClusters in each layer and
0049   // among other the information to compute the SimCluster-To-LayerCluster score. It is one of the two objects that
0050   // build the output of the makeConnections function.
0051   // lcsInSimCluster[scId][layerId]
0052   ticl::simClusterToLayerCluster lcsInSimCluster;
0053   lcsInSimCluster.resize(nSimClusters);
0054   for (unsigned int i = 0; i < nSimClusters; ++i) {
0055     lcsInSimCluster[i].resize(layers_ * 2);
0056     for (unsigned int j = 0; j < layers_ * 2; ++j) {
0057       lcsInSimCluster[i][j].simClusterId = i;
0058       lcsInSimCluster[i][j].energy = 0.f;
0059       lcsInSimCluster[i][j].hits_and_fractions.clear();
0060     }
0061   }
0062
0063   // Fill detIdToSimClusterId_Map and update lcsInSimCluster
0064   // The detIdToSimClusterId_Map is used to connect a hit Detid (key) with all the SimClusters that
0065   // contributed to that hit by storing the SimCluster id and the fraction of the hit. Observe here
0066   // that in contrast to the CaloParticle case there is no merging and summing of the fractions, which
0067   // in the CaloParticle's case was necessary due to the multiple SimClusters of a single CaloParticle.
0068   std::unordered_map<DetId, std::vector<ticl::detIdInfoInCluster>> detIdToSimClusterId_Map;
0069   for (const auto& scId : sCIndices) {
0070     std::vector<std::pair<uint32_t, float>> hits_and_fractions = simClusters[scId].hits_and_fractions();
0071     if constexpr (std::is_same_v<HIT, HGCRecHit>)
0072       hits_and_fractions = simClusters[scId].endcap_hits_and_fractions();
0073     else
0074       hits_and_fractions = simClusters[scId].barrel_hits_and_fractions();
0075     for (const auto& it_haf : hits_and_fractions) {
0076       const auto hitid = (it_haf.first);
0077       unsigned int scLayerId = recHitTools_->getLayer(hitid);
0078       if constexpr (std::is_same_v<HIT, HGCRecHit>)
0079         scLayerId += layers_ * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
0080       const auto itcheck = hitMap_->find(hitid);
0081       if (itcheck != hitMap_->end()) {
0082         auto hit_find_it = detIdToSimClusterId_Map.find(hitid);
0083         if (hit_find_it == detIdToSimClusterId_Map.end()) {
0084           detIdToSimClusterId_Map[hitid] = std::vector<ticl::detIdInfoInCluster>();
0085         }
0086         detIdToSimClusterId_Map[hitid].emplace_back(scId, it_haf.second);
0087         const HIT* hit = hits_[itcheck->second];
0088         lcsInSimCluster[scId][scLayerId].energy += it_haf.second * hit->energy();
0089         lcsInSimCluster[scId][scLayerId].hits_and_fractions.emplace_back(hitid, it_haf.second);
0090       }
0091     }
0092   }  // end of loop over SimClusters
0093
0094 #ifdef EDM_ML_DEBUG
0095   LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0096       << "lcsInSimCluster INFO (Only SimCluster filled at the moment)" << std::endl;
0097   LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "    # of clusters :          " << nLayerClusters << std::endl;
0098   for (size_t sc = 0; sc < lcsInSimCluster.size(); ++sc) {
0099     LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "For SimCluster Idx: " << sc << " we have: " << std::endl;
0100     for (size_t sclay = 0; sclay < lcsInSimCluster[sc].size(); ++sclay) {
0101       LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "  On Layer: " << sclay << " we have:" << std::endl;
0102       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0103           << "    SimClusterIdx: " << lcsInSimCluster[sc][sclay].simClusterId << std::endl;
0104       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0105           << "    Energy:          " << lcsInSimCluster[sc][sclay].energy << std::endl;
0106       double tot_energy = 0.;
0107       for (auto const& haf : lcsInSimCluster[sc][sclay].hits_and_fractions) {
0108         const HIT* hit = hits_[hitMap_->at(haf.first)];
0109         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "      Hits/fraction/energy: " << (uint32_t)haf.first << "/"
0110                                                       << haf.second << "/" << haf.second * hit->energy() << std::endl;
0111         tot_energy += haf.second * hit->energy();
0112       }
0113       LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "    Tot Sum haf: " << tot_energy << std::endl;
0114       for (auto const& lc : lcsInSimCluster[sc][sclay].layerClusterIdToEnergyAndScore) {
0115         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "      lcIdx/energy/score: " << lc.first << "/"
0116                                                       << lc.second.first << "/" << lc.second.second << std::endl;
0117       }
0118     }
0119   }
0120
0121   LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "detIdToSimClusterId_Map INFO" << std::endl;
0122   for (auto const& sc : detIdToSimClusterId_Map) {
0123     LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0124         << "For detId: " << (uint32_t)sc.first
0125         << " we have found the following connections with SimClusters:" << std::endl;
0126     // At this point here if you activate the printing you will notice cases where in a
0127     // specific detId there are more that one SimClusters contributing with fractions less than 1.
0128     // This is important since it effects the score computation, since the fraction is also in the
0129     // denominator of the score formula.
0130     const HIT* hit = hits_[hitMap_->at(sc.first)];
0131     for (auto const& sclu : sc.second) {
0132       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0133           << "  SimCluster Id: " << sclu.clusterId << " with fraction: " << sclu.fraction
0134           << " and energy: " << sclu.fraction * hit->energy() << std::endl;
0135     }
0136   }
0137 #endif
0138
0139   // Fill detIdToLayerClusterId_Map and scsInLayerCluster; update lcsInSimCluster
0140   // The detIdToLayerClusterId_Map is used to connect a hit Detid (key) with all the LayerClusters that
0141   // contributed to that hit by storing the LayerCluster id and the fraction of the corresponding hit.
0142   std::unordered_map<DetId, std::vector<ticl::detIdInfoInCluster>> detIdToLayerClusterId_Map;
0143   // scsInLayerCluster together with lcsInSimCluster are the two objects that are used to build the
0144   // output of the makeConnections function. scsInLayerCluster connects a LayerCluster with
0145   // all the SimClusters that share at least one cell with the LayerCluster and for each pair (LC,SC)
0146   // it stores the score.
0147   ticl::layerClusterToSimCluster scsInLayerCluster;  //[lcId][scId]->(score)
0148   scsInLayerCluster.resize(nLayerClusters);
0149
0150   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0151     const std::vector<std::pair<DetId, float>>& hits_and_fractions = clusters[lcId].hitsAndFractions();
0152     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0153     const auto firstHitDetId = hits_and_fractions[0].first;
0154     int lcLayerId = recHitTools_->getLayer(firstHitDetId);
0155     if constexpr (std::is_same_v<HIT, HGCRecHit>)
0156       lcLayerId += layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
0157     for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0158       const auto rh_detid = hits_and_fractions[hitId].first;
0159       const auto rhFraction = hits_and_fractions[hitId].second;
0160
0161       auto hit_find_in_LC = detIdToLayerClusterId_Map.find(rh_detid);
0162       if (hit_find_in_LC == detIdToLayerClusterId_Map.end()) {
0163         detIdToLayerClusterId_Map[rh_detid] = std::vector<ticl::detIdInfoInCluster>();
0164       }
0165       detIdToLayerClusterId_Map[rh_detid].emplace_back(lcId, rhFraction);
0166
0167       auto hit_find_in_SC = detIdToSimClusterId_Map.find(rh_detid);
0168
0169       if (hit_find_in_SC != detIdToSimClusterId_Map.end()) {
0170         const auto itcheck = hitMap_->find(rh_detid);
0171         const HIT* hit = hits_[itcheck->second];
0172         //Loops through all the simclusters that have the layer cluster rechit under study
0173         //Here is time to update the lcsInSimCluster and connect the SimCluster with all
0174         //the layer clusters that have the current rechit detid matched.
0175         for (auto& h : hit_find_in_SC->second) {
0176           //lcsInSimCluster[simclusterId][layerId][layerclusterId]-> (energy,score)
0177           //SC_i - > LC_j, LC_k, ...
0178           lcsInSimCluster[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[lcId].first +=
0179               h.fraction * hit->energy();
0180           //LC_i -> SC_j, SC_k, ...
0181           scsInLayerCluster[lcId].emplace_back(h.clusterId, 0.f);
0182         }
0183       }
0184     }  // End loop over hits on a LayerCluster
0185   }  // End of loop over LayerClusters
0186
0187 #ifdef EDM_ML_DEBUG
0188   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0189     const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
0190     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0191     const auto firstHitDetId = hits_and_fractions[0].first;
0192     int lcLayerId = recHitTools_->getLayer(firstHitDetId);
0193     if constexpr (std::is_same_v<HIT, HGCRecHit>)
0194       lcLayerId += layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
0195     // This vector will store, for each hit in the Layercluster, the index of
0196     // the SimCluster that contributed the most, in terms of energy, to it.
0197     // Special values are:
0198     //
0199     // -2  --> the reconstruction fraction of the RecHit is 0 (used in the past to monitor Halo Hits)
0200     // -3  --> same as before with the added condition that no SimCluster has been linked to this RecHit
0201     // -1  --> the reco fraction is >0, but no SimCluster has been linked to it
0202     // >=0 --> index of the linked SimCluster
0203     std::vector<int> hitsToSimClusterId(numberOfHitsInLC);
0204     // This will store the index of the SimCluster linked to the LayerCluster that has the most number of hits in common.
0205     int maxSCId_byNumberOfHits = -1;
0206     // This will store the maximum number of shared hits between a Layercluster and a SimCluster
0207     unsigned int maxSCNumberOfHitsInLC = 0;
0208     // This will store the index of the SimCluster linked to the LayerCluster that has the most energy in common.
0209     int maxSCId_byEnergy = -1;
0210     // This will store the maximum number of shared energy between a Layercluster and a SimCluster
0211     float maxEnergySharedLCandSC = 0.f;
0212     // This will store the fraction of the LayerCluster energy shared with the best(energy) SimCluster: e_shared/lc_energy
0213     float energyFractionOfLCinSC = 0.f;
0214     // This will store the fraction of the SimCluster energy shared with the LayerCluster: e_shared/sc_energy
0215     float energyFractionOfSCinLC = 0.f;
0216     std::unordered_map<unsigned, unsigned> occurrencesSCinLC;
0217     unsigned int numberOfNoiseHitsInLC = 0;
0218     std::unordered_map<unsigned, float> SCEnergyInLC;
0219
0220     for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0221       const auto rh_detid = hits_and_fractions[hitId].first;
0222       const auto rhFraction = hits_and_fractions[hitId].second;
0223
0224       auto hit_find_in_SC = detIdToSimClusterId_Map.find(rh_detid);
0225
0226       // if the fraction is zero or the hit does not belong to any simcluster
0227       // set the caloparticleId for the hit to -1 this will
0228       // contribute to the number of noise hits
0229
0230       // MR Remove the case in which the fraction is 0, since this could be a
0231       // real hit that has been marked as halo.
0232       if (rhFraction == 0.) {
0233         hitsToSimClusterId[hitId] = -2;
0234       }
0235       //Now check if there are SimClusters linked to this rechit of the layercluster
0236       if (hit_find_in_SC == detIdToSimClusterId_Map.end()) {
0237         hitsToSimClusterId[hitId] -= 1;
0238       } else {
0239         const auto itcheck = hitMap_->find(rh_detid);
0240         const HIT* hit = hits_[itcheck->second];
0241         auto maxSCEnergyInLC = 0.f;
0242         auto maxSCId = -1;
0243         //Loop through all the linked SimClusters
0244         for (auto& h : hit_find_in_SC->second) {
0245           SCEnergyInLC[h.clusterId] += h.fraction * hit->energy();
0246           // Keep track of which SimCluster ccontributed the most, in terms
0247           // of energy, to this specific LayerCluster.
0248           if (SCEnergyInLC[h.clusterId] > maxSCEnergyInLC) {
0249             maxSCEnergyInLC = SCEnergyInLC[h.clusterId];
0250             maxSCId = h.clusterId;
0251           }
0252         }
0253         hitsToSimClusterId[hitId] = maxSCId;
0254       }
0255     }  // End loop over hits on a LayerCluster
0256
0257     for (const auto& c : hitsToSimClusterId) {
0258       if (c < 0) {
0259         numberOfNoiseHitsInLC++;
0260       } else {
0261         occurrencesSCinLC[c]++;
0262       }
0263     }
0264
0265     for (const auto& c : occurrencesSCinLC) {
0266       if (c.second > maxSCNumberOfHitsInLC) {
0267         maxSCId_byNumberOfHits = c.first;
0268         maxSCNumberOfHitsInLC = c.second;
0269       }
0270     }
0271
0272     for (const auto& c : SCEnergyInLC) {
0273       if (c.second > maxEnergySharedLCandSC) {
0274         maxSCId_byEnergy = c.first;
0275         maxEnergySharedLCandSC = c.second;
0276       }
0277     }
0278
0279     float totalSCEnergyOnLayer = 0.f;
0280     if (maxSCId_byEnergy >= 0) {
0281       totalSCEnergyOnLayer = lcsInSimCluster[maxSCId_byEnergy][lcLayerId].energy;
0282       energyFractionOfSCinLC = maxEnergySharedLCandSC / totalSCEnergyOnLayer;
0283       if (clusters[lcId].energy() > 0.f) {
0284         energyFractionOfLCinSC = maxEnergySharedLCandSC / clusters[lcId].energy();
0285       }
0286     }
0287
0288     LogDebug("LCToSCAssociatorByEnergyScoreImpl") << std::setw(10) << "LayerId:"
0289                                                   << "\t" << std::setw(12) << "layerCluster"
0290                                                   << "\t" << std::setw(10) << "lc energy"
0291                                                   << "\t" << std::setw(5) << "nhits"
0292                                                   << "\t" << std::setw(12) << "noise hits"
0293                                                   << "\t" << std::setw(22) << "maxSCId_byNumberOfHits"
0294                                                   << "\t" << std::setw(8) << "nhitsSC"
0295                                                   << "\t" << std::setw(13) << "maxSCId_byEnergy"
0296                                                   << "\t" << std::setw(20) << "maxEnergySharedLCandSC"
0297                                                   << "\t" << std::setw(22) << "totalSCEnergyOnLayer"
0298                                                   << "\t" << std::setw(22) << "energyFractionOfLCinSC"
0299                                                   << "\t" << std::setw(25) << "energyFractionOfSCinLC"
0300                                                   << "\t"
0301                                                   << "\n";
0302     LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0303         << std::setw(10) << lcLayerId << "\t" << std::setw(12) << lcId << "\t" << std::setw(10)
0304         << clusters[lcId].energy() << "\t" << std::setw(5) << numberOfHitsInLC << "\t" << std::setw(12)
0305         << numberOfNoiseHitsInLC << "\t" << std::setw(22) << maxSCId_byNumberOfHits << "\t" << std::setw(8)
0306         << maxSCNumberOfHitsInLC << "\t" << std::setw(13) << maxSCId_byEnergy << "\t" << std::setw(20)
0307         << maxEnergySharedLCandSC << "\t" << std::setw(22) << totalSCEnergyOnLayer << "\t" << std::setw(22)
0308         << energyFractionOfLCinSC << "\t" << std::setw(25) << energyFractionOfSCinLC << "\n";
0309   }  // End of loop over LayerClusters
0310
0311   LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0312       << "Improved lcsInSimCluster INFO (Now containing the linked layer clusters id and energy - score still empty)"
0313       << std::endl;
0314   for (size_t sc = 0; sc < lcsInSimCluster.size(); ++sc) {
0315     LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "For SimCluster Idx: " << sc << " we have: " << std::endl;
0316     for (size_t sclay = 0; sclay < lcsInSimCluster[sc].size(); ++sclay) {
0317       LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "  On Layer: " << sclay << " we have:" << std::endl;
0318       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0319           << "    SimClusterIdx: " << lcsInSimCluster[sc][sclay].simClusterId << std::endl;
0320       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0321           << "    Energy:          " << lcsInSimCluster[sc][sclay].energy << std::endl;
0322       double tot_energy = 0.;
0323       for (auto const& haf : lcsInSimCluster[sc][sclay].hits_and_fractions) {
0324         const HIT* hit = hits_[hitMap_->at(haf.first)];
0325         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "      Hits/fraction/energy: " << (uint32_t)haf.first << "/"
0326                                                       << haf.second << "/" << haf.second * hit->energy() << std::endl;
0327         tot_energy += haf.second * hit->energy();
0328       }
0329       LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "    Tot Sum haf: " << tot_energy << std::endl;
0330       for (auto const& lc : lcsInSimCluster[sc][sclay].layerClusterIdToEnergyAndScore) {
0331         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "      lcIdx/energy/score: " << lc.first << "/"
0332                                                       << lc.second.first << "/" << lc.second.second << std::endl;
0333       }
0334     }
0335   }
0336
0337   LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "Improved detIdToSimClusterId_Map INFO" << std::endl;
0338   for (auto const& sc : detIdToSimClusterId_Map) {
0339     const HIT* hit = hits_[hitMap_->at(sc.first)];
0340     LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0341         << "For detId: " << (uint32_t)sc.first
0342         << " we have found the following connections with SimClusters:" << std::endl;
0343     for (auto const& sclu : sc.second) {
0344       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0345           << "  SimCluster Id: " << sclu.clusterId << " with fraction: " << sclu.fraction
0346           << " and energy: " << sclu.fraction * hit->energy() << std::endl;
0347     }
0348   }
0349 #endif
0350
0351   // Update scsInLayerCluster; compute the score LayerCluster-to-SimCluster,
0352   // together with the returned AssociationMap
0353   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0354     // The simclusters contributing to the layer clusters should already be unique.
0355     // find the unique simclusters id contributing to the layer clusters
0356     std::sort(scsInLayerCluster[lcId].begin(), scsInLayerCluster[lcId].end());
0357     auto last = std::unique(scsInLayerCluster[lcId].begin(), scsInLayerCluster[lcId].end());
0358     scsInLayerCluster[lcId].erase(last, scsInLayerCluster[lcId].end());
0359     const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
0360     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0361     // If a reconstructed LayerCluster has energy 0 but is linked to a
0362     // SimCluster, assigned score 1
0363     if (clusters[lcId].energy() == 0. && !scsInLayerCluster[lcId].empty()) {
0364       for (auto& scPair : scsInLayerCluster[lcId]) {
0365         scPair.second = 1.;
0366         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "layerClusterId : \t " << lcId << "\t SC id : \t"
0367                                                       << scPair.first << "\t score \t " << scPair.second << "\n";
0368       }
0369       continue;
0370     }
0371
0372     // Compute the correct normalization.
0373     // It is the inverse of the denominator of the LCToSC score formula. Observe that this is the sum of the squares.
0374     float invLayerClusterEnergyWeight = 0.f;
0375     for (auto const& haf : hits_and_fractions) {
0376       const HIT* hit = hits_[hitMap_->at(haf.first)];
0377       invLayerClusterEnergyWeight += (haf.second * hit->energy()) * (haf.second * hit->energy());
0378     }
0379     invLayerClusterEnergyWeight = 1.f / invLayerClusterEnergyWeight;
0380     for (unsigned int i = 0; i < numberOfHitsInLC; ++i) {
0381       DetId rh_detid = hits_and_fractions[i].first;
0382       float rhFraction = hits_and_fractions[i].second;
0383
0384       bool hitWithSC = (detIdToSimClusterId_Map.find(rh_detid) != detIdToSimClusterId_Map.end());
0385
0386       auto itcheck = hitMap_->find(rh_detid);
0387       const HIT* hit = hits_[itcheck->second];
0388       float hitEnergyWeight = hit->energy() * hit->energy();
0389
0390       for (auto& scPair : scsInLayerCluster[lcId]) {
0391         float scFraction = 0.f;
0392         if (hitWithSC) {
0393           auto findHitIt = std::find(detIdToSimClusterId_Map[rh_detid].begin(),
0394                                      detIdToSimClusterId_Map[rh_detid].end(),
0395                                      ticl::detIdInfoInCluster{scPair.first, 0.f});
0396           if (findHitIt != detIdToSimClusterId_Map[rh_detid].end())
0397             scFraction = findHitIt->fraction;
0398         }
0399         scPair.second += std::min(std::pow(rhFraction - scFraction, 2), std::pow(rhFraction, 2)) * hitEnergyWeight *
0400                          invLayerClusterEnergyWeight;
0401 #ifdef EDM_ML_DEBUG
0402         LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0403             << "rh_detid:\t" << (uint32_t)rh_detid << "\tlayerClusterId:\t" << lcId << "\t"
0404             << "rhfraction,scfraction:\t" << rhFraction << ", " << scFraction << "\t"
0405             << "hitEnergyWeight:\t" << hitEnergyWeight << "\t"
0406             << "current score:\t" << scPair.second << "\t"
0407             << "invLayerClusterEnergyWeight:\t" << invLayerClusterEnergyWeight << "\n";
0408 #endif
0409       }
0410     }  // End of loop over Hits within a LayerCluster
0411 #ifdef EDM_ML_DEBUG
0412     if (scsInLayerCluster[lcId].empty())
0413       LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "layerCluster Id: \t" << lcId << "\tSC id:\t-1 "
0414                                                     << "\t score \t-1"
0415                                                     << "\n";
0416 #endif
0417   }  // End of loop over LayerClusters
0418
0419   // Compute the SimCluster-To-LayerCluster score
0420   for (const auto& scId : sCIndices) {
0421     for (unsigned int layerId = 0; layerId < layers_ * 2; ++layerId) {
0422       unsigned int SCNumberOfHits = lcsInSimCluster[scId][layerId].hits_and_fractions.size();
0423       if (SCNumberOfHits == 0)
0424         continue;
0425 #ifdef EDM_ML_DEBUG
0426       int lcWithMaxEnergyInSC = -1;
0427       //energy of the most energetic LC from all that were linked to SC
0428       float maxEnergyLCinSC = 0.f;
0429       //Energy of the SC scId on layer layerId that was reconstructed.
0430       float SCenergy = lcsInSimCluster[scId][layerId].energy;
0431       //most energetic LC from all LCs linked to SC over SC energy.
0432       float SCEnergyFractionInLC = 0.f;
0433       for (auto& lc : lcsInSimCluster[scId][layerId].layerClusterIdToEnergyAndScore) {
0434         if (lc.second.first > maxEnergyLCinSC) {
0435           maxEnergyLCinSC = lc.second.first;
0436           lcWithMaxEnergyInSC = lc.first;
0437         }
0438       }
0439       if (SCenergy > 0.f)
0440         SCEnergyFractionInLC = maxEnergyLCinSC / SCenergy;
0441
0442       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0443           << std::setw(8) << "LayerId:\t" << std::setw(12) << "simcluster\t" << std::setw(15) << "sc total energy\t"
0444           << std::setw(15) << "scEnergyOnLayer\t" << std::setw(14) << "SCNhitsOnLayer\t" << std::setw(18)
0445           << "lcWithMaxEnergyInSC\t" << std::setw(15) << "maxEnergyLCinSC\t" << std::setw(20) << "SCEnergyFractionInLC"
0446           << "\n";
0447       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0448           << std::setw(8) << layerId << "\t" << std::setw(12) << scId << "\t" << std::setw(15)
0449           << simClusters[scId].energy() << "\t" << std::setw(15) << SCenergy << "\t" << std::setw(14) << SCNumberOfHits
0450           << "\t" << std::setw(18) << lcWithMaxEnergyInSC << "\t" << std::setw(15) << maxEnergyLCinSC << "\t"
0451           << std::setw(20) << SCEnergyFractionInLC << "\n";
0452 #endif
0453       // Compute the correct normalization. Observe that this is the sum of the squares.
0454       float invSCEnergyWeight = 0.f;
0455       for (auto const& haf : lcsInSimCluster[scId][layerId].hits_and_fractions) {
0456         const HIT* hit = hits_[hitMap_->at(haf.first)];
0457         invSCEnergyWeight += std::pow(haf.second * hit->energy(), 2);
0458       }
0459       invSCEnergyWeight = 1.f / invSCEnergyWeight;
0460       for (unsigned int i = 0; i < SCNumberOfHits; ++i) {
0461         auto& sc_hitDetId = lcsInSimCluster[scId][layerId].hits_and_fractions[i].first;
0462         auto& scFraction = lcsInSimCluster[scId][layerId].hits_and_fractions[i].second;
0463
0464         bool hitWithLC = false;
0465         if (scFraction == 0.f)
0466           continue;  //hopefully this should never happen
0467         auto hit_find_in_LC = detIdToLayerClusterId_Map.find(sc_hitDetId);
0468         if (hit_find_in_LC != detIdToLayerClusterId_Map.end())
0469           hitWithLC = true;
0470         auto itcheck = hitMap_->find(sc_hitDetId);
0471         const HIT* hit = hits_[itcheck->second];
0472         float hitEnergyWeight = hit->energy() * hit->energy();
0473         for (auto& lcPair : lcsInSimCluster[scId][layerId].layerClusterIdToEnergyAndScore) {
0474           unsigned int layerClusterId = lcPair.first;
0475           float lcFraction = 0.f;
0476
0477           if (hitWithLC) {
0478             auto findHitIt = std::find(detIdToLayerClusterId_Map[sc_hitDetId].begin(),
0479                                        detIdToLayerClusterId_Map[sc_hitDetId].end(),
0480                                        ticl::detIdInfoInCluster{layerClusterId, 0.f});
0481             if (findHitIt != detIdToLayerClusterId_Map[sc_hitDetId].end())
0482               lcFraction = findHitIt->fraction;
0483           }
0484           lcPair.second.second += std::min(std::pow(lcFraction - scFraction, 2), std::pow(scFraction, 2)) *
0485                                   hitEnergyWeight * invSCEnergyWeight;
0486 #ifdef EDM_ML_DEBUG
0487           LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0488               << "scDetId:\t" << (uint32_t)sc_hitDetId << "\tlayerClusterId:\t" << layerClusterId << "\t"
0489               << "lcfraction,scfraction:\t" << lcFraction << ", " << scFraction << "\t"
0490               << "hitEnergyWeight:\t" << hitEnergyWeight << "\t"
0491               << "current score:\t" << lcPair.second.second << "\t"
0492               << "invSCEnergyWeight:\t" << invSCEnergyWeight << "\n";
0493 #endif
0494         }  // End of loop over LayerClusters linked to hits of this SimCluster
0495       }  // End of loop over hits of SimCluster on a Layer
0496 #ifdef EDM_ML_DEBUG
0497       if (lcsInSimCluster[scId][layerId].layerClusterIdToEnergyAndScore.empty())
0498         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "SC Id: \t" << scId << "\tLC id:\t-1 "
0499                                                       << "\t score \t-1"
0500                                                       << "\n";
0501
0502       for (const auto& lcPair : lcsInSimCluster[scId][layerId].layerClusterIdToEnergyAndScore) {
0503         LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0504             << "SC Id: \t" << scId << "\t LC id: \t" << lcPair.first << "\t score \t" << lcPair.second.second
0505             << "\t shared energy:\t" << lcPair.second.first << "\t shared energy fraction:\t"
0506             << (lcPair.second.first / SCenergy) << "\n";
0507       }
0508 #endif
0509     }  // End of loop over layers
0510   }  // End of loop over SimClusters
0511
0512   return {scsInLayerCluster, lcsInSimCluster};
0513 }
0514
0515 template <typename HIT>
0516 ticl::RecoToSimCollectionWithSimClusters LCToSCAssociatorByEnergyScoreImpl<HIT>::associateRecoToSim(
0517     const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<SimClusterCollection>& sCCH) const {
0518   ticl::RecoToSimCollectionWithSimClusters returnValue(productGetter_);
0519   const auto& links = makeConnections(cCCH, sCCH);
0520
0521   const auto& scsInLayerCluster = std::get<0>(links);
0522   for (size_t lcId = 0; lcId < scsInLayerCluster.size(); ++lcId) {
0523     for (auto& scPair : scsInLayerCluster[lcId]) {
0524       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0525           << "layerCluster Id: \t" << lcId << "\t SC id: \t" << scPair.first << "\t score \t" << scPair.second << "\n";
0526       // Fill AssociationMap
0527       returnValue.insert(edm::Ref<reco::CaloClusterCollection>(cCCH, lcId),  // Ref to LC
0528                          std::make_pair(edm::Ref<SimClusterCollection>(sCCH, scPair.first),
0529                                         scPair.second)  // Pair <Ref to SC, score>
0530       );
0531     }
0532   }
0533   return returnValue;
0534 }
0535
0536 template <typename HIT>
0537 ticl::SimToRecoCollectionWithSimClusters LCToSCAssociatorByEnergyScoreImpl<HIT>::associateSimToReco(
0538     const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<SimClusterCollection>& sCCH) const {
0539   ticl::SimToRecoCollectionWithSimClusters returnValue(productGetter_);
0540   const auto& links = makeConnections(cCCH, sCCH);
0541   const auto& lcsInSimCluster = std::get<1>(links);
0542   for (size_t scId = 0; scId < lcsInSimCluster.size(); ++scId) {
0543     for (size_t layerId = 0; layerId < lcsInSimCluster[scId].size(); ++layerId) {
0544       for (auto& lcPair : lcsInSimCluster[scId][layerId].layerClusterIdToEnergyAndScore) {
0545         returnValue.insert(
0546             edm::Ref<SimClusterCollection>(sCCH, scId),                                // Ref to SC
0547             std::make_pair(edm::Ref<reco::CaloClusterCollection>(cCCH, lcPair.first),  // Pair <Ref to LC,
0548                            std::make_pair(lcPair.second.first, lcPair.second.second))  // pair <energy, score> >
0549         );
0550       }
0551     }
0552   }
0553   return returnValue;
0554 }
0555
0556 template class LCToSCAssociatorByEnergyScoreImpl<HGCRecHit>;
0557 template class LCToSCAssociatorByEnergyScoreImpl<reco::PFRecHit>;