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