Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-27 23:06:24

0001 // Original Author: Leonardo Cristella
0002 //
0003 
0004 #include "MultiClusterAssociatorByEnergyScoreImpl.h"
0005 
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0008 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0009 #include "SimCalorimetry/HGCalAssociatorProducers/interface/AssociatorTools.h"
0010 
0011 #include <cfloat>
0012 
0013 MultiClusterAssociatorByEnergyScoreImpl::MultiClusterAssociatorByEnergyScoreImpl(
0014     edm::EDProductGetter const& productGetter,
0015     bool hardScatterOnly,
0016     std::shared_ptr<hgcal::RecHitTools> recHitTools,
0017     const std::unordered_map<DetId, const HGCRecHit*>*& hitMap)
0018     : hardScatterOnly_(hardScatterOnly), recHitTools_(recHitTools), hitMap_(hitMap), productGetter_(&productGetter) {
0019   layers_ = recHitTools_->lastLayerBH();
0020 }
0021 
0022 hgcal::association MultiClusterAssociatorByEnergyScoreImpl::makeConnections(
0023     const edm::Handle<reco::HGCalMultiClusterCollection>& mCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0024   // 1. Extract collections and filter CaloParticles, if required
0025   const auto& mClusters = *mCCH.product();
0026   const auto& caloParticles = *cPCH.product();
0027   auto nMultiClusters = mClusters.size();
0028   //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution.
0029   std::vector<size_t> cPIndices;
0030   //Consider CaloParticles coming from the hard scatterer
0031   //excluding the PU contribution and save the indices.
0032   removeCPFromPU(caloParticles, cPIndices, false);
0033   auto nCaloParticles = cPIndices.size();
0034 
0035   std::vector<size_t> cPSelectedIndices;
0036   removeCPFromPU(caloParticles, cPSelectedIndices, true);
0037 
0038   //cPOnLayer[caloparticle][layer]
0039   //This defines a "caloParticle on layer" concept. It is only filled in case
0040   //that caloParticle has a reconstructed hit related via detid. So, a cPOnLayer[i][j] connects a
0041   //specific caloParticle i in layer j with:
0042   //1. the sum of all recHits energy times fraction of the relevant simHit in layer j related to that caloParticle i.
0043   //2. the hits and fractions of that caloParticle i in layer j.
0044   //3. the layer clusters with matched recHit id.
0045   hgcal::caloParticleToMultiCluster cPOnLayer;
0046   cPOnLayer.resize(nCaloParticles);
0047   for (unsigned int i = 0; i < nCaloParticles; ++i) {
0048     auto cpIndex = cPIndices[i];
0049     cPOnLayer[cpIndex].resize(layers_ * 2);
0050     for (unsigned int j = 0; j < layers_ * 2; ++j) {
0051       cPOnLayer[cpIndex][j].caloParticleId = cpIndex;
0052       cPOnLayer[cpIndex][j].energy = 0.f;
0053       cPOnLayer[cpIndex][j].hits_and_fractions.clear();
0054     }
0055   }
0056 
0057   std::unordered_map<DetId, std::vector<hgcal::detIdInfoInCluster>> detIdToCaloParticleId_Map;
0058   // Fill detIdToCaloParticleId_Map and update cPOnLayer
0059   for (const auto& cpId : cPIndices) {
0060     //take sim clusters
0061     const SimClusterRefVector& simClusterRefVector = caloParticles[cpId].simClusters();
0062     //loop through sim clusters
0063     for (const auto& it_sc : simClusterRefVector) {
0064       const SimCluster& simCluster = (*(it_sc));
0065       const auto& hits_and_fractions = simCluster.hits_and_fractions();
0066       for (const auto& it_haf : hits_and_fractions) {
0067         const auto hitid = (it_haf.first);
0068         const auto cpLayerId =
0069             recHitTools_->getLayerWithOffset(hitid) + layers_ * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
0070         const auto itcheck = hitMap_->find(hitid);
0071         if (itcheck != hitMap_->end()) {
0072           //Since the current hit from sim cluster has a reconstructed hit with the same detid,
0073           //make a map that will connect a detid with:
0074           //1. the caloParticles that have a simcluster with sim hits in that cell via caloParticle id.
0075           //2. the sum of all simHits fractions that contributes to that detid.
0076           //So, keep in mind that in case of multiple caloParticles contributing in the same cell
0077           //the fraction is the sum over all caloParticles. So, something like:
0078           //detid: (caloParticle 1, sum of hits fractions in that detid over all cp) , (caloParticle 2, sum of hits fractions in that detid over all cp), (caloParticle 3, sum of hits fractions in that detid over all cp) ...
0079           auto hit_find_it = detIdToCaloParticleId_Map.find(hitid);
0080           if (hit_find_it == detIdToCaloParticleId_Map.end()) {
0081             detIdToCaloParticleId_Map[hitid] = std::vector<hgcal::detIdInfoInCluster>();
0082             detIdToCaloParticleId_Map[hitid].emplace_back(cpId, it_haf.second);
0083           } else {
0084             auto findHitIt = std::find(detIdToCaloParticleId_Map[hitid].begin(),
0085                                        detIdToCaloParticleId_Map[hitid].end(),
0086                                        hgcal::detIdInfoInCluster{cpId, it_haf.second});
0087             if (findHitIt != detIdToCaloParticleId_Map[hitid].end()) {
0088               findHitIt->fraction += it_haf.second;
0089             } else {
0090               detIdToCaloParticleId_Map[hitid].emplace_back(cpId, it_haf.second);
0091             }
0092           }
0093           const HGCRecHit* hit = itcheck->second;
0094           //Since the current hit from sim cluster has a reconstructed hit with the same detid,
0095           //fill the cPOnLayer[caloparticle][layer] object with energy (sum of all recHits energy times fraction
0096           //of the relevant simHit) and keep the hit (detid and fraction) that contributed.
0097           cPOnLayer[cpId][cpLayerId].energy += it_haf.second * hit->energy();
0098           // We need to compress the hits and fractions in order to have a
0099           // reasonable score between CP and LC. Imagine, for example, that a
0100           // CP has detID X used by 2 SimClusters with different fractions. If
0101           // a single LC uses X with fraction 1 and is compared to the 2
0102           // contributions separately, it will be assigned a score != 0, which
0103           // is wrong.
0104           auto& haf = cPOnLayer[cpId][cpLayerId].hits_and_fractions;
0105           auto found = std::find_if(
0106               std::begin(haf), std::end(haf), [&hitid](const std::pair<DetId, float>& v) { return v.first == hitid; });
0107           if (found != haf.end()) {
0108             found->second += it_haf.second;
0109           } else {
0110             cPOnLayer[cpId][cpLayerId].hits_and_fractions.emplace_back(hitid, it_haf.second);
0111           }
0112         }
0113       }  // end of loop through simHits
0114     }    // end of loop through simclusters
0115   }      // end of loop through caloParticles
0116 
0117 #ifdef EDM_ML_DEBUG
0118   LogDebug("MultiClusterAssociatorByEnergyScoreImpl") << "cPOnLayer INFO" << std::endl;
0119   for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) {
0120     LogDebug("MultiClusterAssociatorByEnergyScoreImpl") << "For CaloParticle Idx: " << cp << " we have: " << std::endl;
0121     for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) {
0122       LogDebug("MultiClusterAssociatorByEnergyScoreImpl") << "  On Layer: " << cpp << " we have:" << std::endl;
0123       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0124           << "    CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl;
0125       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0126           << "    Energy:          " << cPOnLayer[cp][cpp].energy << std::endl;
0127       double tot_energy = 0.;
0128       for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) {
0129         LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0130             << "      Hits/fraction/energy: " << (uint32_t)haf.first << "/" << haf.second << "/"
0131             << haf.second * hitMap_->at(haf.first)->energy() << std::endl;
0132         tot_energy += haf.second * hitMap_->at(haf.first)->energy();
0133       }
0134       LogDebug("MultiClusterAssociatorByEnergyScoreImpl") << "    Tot Sum haf: " << tot_energy << std::endl;
0135       for (auto const& mc : cPOnLayer[cp][cpp].multiClusterIdToEnergyAndScore) {
0136         LogDebug("MultiClusterAssociatorByEnergyScoreImpl") << "      mcIdx/energy/score: " << mc.first << "/"
0137                                                             << mc.second.first << "/" << mc.second.second << std::endl;
0138       }
0139     }
0140   }
0141 
0142   LogDebug("MultiClusterAssociatorByEnergyScoreImpl") << "detIdToCaloParticleId_Map INFO" << std::endl;
0143   for (auto const& cp : detIdToCaloParticleId_Map) {
0144     LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0145         << "For detId: " << (uint32_t)cp.first
0146         << " we have found the following connections with CaloParticles:" << std::endl;
0147     for (auto const& cpp : cp.second) {
0148       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0149           << "  CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction
0150           << " and energy: " << cpp.fraction * hitMap_->at(cp.first)->energy() << std::endl;
0151     }
0152   }
0153 #endif
0154 
0155   // Fill detIdToMultiClusterId_Map and cpsInMultiCluster; update cPOnLayer
0156   std::unordered_map<DetId, std::vector<hgcal::detIdInfoInMultiCluster>> detIdToMultiClusterId_Map;
0157 
0158   // this contains the ids of the caloParticles contributing with at least one hit to the multiCluster and the reconstruction error
0159   //cpsInMultiCluster[multicluster][CPids]
0160   //Connects a multiCluster with all related caloParticles.
0161   hgcal::multiClusterToCaloParticle cpsInMultiCluster;
0162   cpsInMultiCluster.resize(nMultiClusters);
0163 
0164   //Loop through multiClusters
0165   for (unsigned int mcId = 0; mcId < nMultiClusters; ++mcId) {
0166     const auto& hits_and_fractions = mClusters[mcId].hitsAndFractions();
0167     if (!hits_and_fractions.empty()) {
0168       std::unordered_map<unsigned, float> CPEnergyInMCL;
0169       int maxCPId_byNumberOfHits = -1;
0170       unsigned int maxCPNumberOfHitsInMCL = 0;
0171       int maxCPId_byEnergy = -1;
0172       float maxEnergySharedMCLandCP = 0.f;
0173       float energyFractionOfMCLinCP = 0.f;
0174       float energyFractionOfCPinMCL = 0.f;
0175 
0176       //In case of matched rechit-simhit, so matched
0177       //caloparticle-layercluster-multicluster, we count and save the number of
0178       //recHits related to the maximum energy CaloParticle out of all
0179       //CaloParticles related to that layer cluster and multiCluster.
0180 
0181       std::unordered_map<unsigned, unsigned> occurrencesCPinMCL;
0182       unsigned int numberOfNoiseHitsInMCL = 0;
0183       unsigned int numberOfHaloHitsInMCL = 0;
0184       unsigned int numberOfHitsInMCL = 0;
0185 
0186       //number of hits related to that cluster
0187       unsigned int numberOfHitsInLC = hits_and_fractions.size();
0188       numberOfHitsInMCL += numberOfHitsInLC;
0189       std::unordered_map<unsigned, float> CPEnergyInLC;
0190 
0191       //hitsToCaloParticleId is a vector of ints, one for each recHit of the
0192       //layer cluster under study. If negative, there is no simHit from any CaloParticle related.
0193       //If positive, at least one CaloParticle has been found with matched simHit.
0194       //In more detail:
0195       // 1. hitsToCaloParticleId[hitId] = -3
0196       //    TN:  These represent Halo Cells(N) that have not been
0197       //    assigned to any CaloParticle (hence the T).
0198       // 2. hitsToCaloParticleId[hitId] = -2
0199       //    FN: There represent Halo Cells(N) that have been assigned
0200       //    to a CaloParticle (hence the F, since those should have not been marked as halo)
0201       // 3. hitsToCaloParticleId[hitId] = -1
0202       //    FP: These represent Real Cells(P) that have not been
0203       //    assigned to any CaloParticle (hence the F, since these are fakes)
0204       // 4. hitsToCaloParticleId[hitId] >= 0
0205       //    TP There represent Real Cells(P) that have been assigned
0206       //    to a CaloParticle (hence the T)
0207 
0208       std::vector<int> hitsToCaloParticleId(numberOfHitsInLC);
0209       //det id of the first hit just to make the lcLayerId variable
0210       //which maps the layers in -z: 0->51 and in +z: 52->103
0211       const auto firstHitDetId = hits_and_fractions[0].first;
0212       int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) +
0213                       layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
0214 
0215       //Loop though the hits of the layer cluster under study
0216       for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0217         const auto rh_detid = hits_and_fractions[hitId].first;
0218         const auto rhFraction = hits_and_fractions[hitId].second;
0219 
0220         //Since the hit is belonging to the layer cluster, it must also be in the recHits map.
0221         const auto itcheck = hitMap_->find(rh_detid);
0222         const auto hit = itcheck->second;
0223 
0224         //Make a map that will connect a detid (that belongs to a recHit of the layer cluster under study,
0225         //no need to save others) with:
0226         //1. the layer clusters that have recHits in that detid
0227         //2. the fraction of the recHit of each layer cluster that contributes to that detid.
0228         //So, something like:
0229         //detid: (layer cluster 1, hit fraction) , (layer cluster 2, hit fraction), (layer cluster 3, hit fraction) ...
0230         //here comparing with the caloParticle map above
0231         auto hit_find_in_LC = detIdToMultiClusterId_Map.find(rh_detid);
0232         if (hit_find_in_LC == detIdToMultiClusterId_Map.end()) {
0233           detIdToMultiClusterId_Map[rh_detid] = std::vector<hgcal::detIdInfoInMultiCluster>();
0234         }
0235         detIdToMultiClusterId_Map[rh_detid].emplace_back(hgcal::detIdInfoInMultiCluster{mcId, mcId, rhFraction});
0236 
0237         // Check whether the recHit of the layer cluster under study has a sim hit in the same cell
0238         auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
0239 
0240         // If the fraction is zero or the hit does not belong to any calo
0241         // particle, set the caloParticleId for the hit to -1 and this will
0242         // contribute to the number of noise hits
0243         if (rhFraction == 0.) {  // this could be a real hit that has been marked as halo
0244           hitsToCaloParticleId[hitId] = -2;
0245           numberOfHaloHitsInMCL++;
0246         }
0247         if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) {
0248           hitsToCaloParticleId[hitId] -= 1;
0249         } else {
0250           auto maxCPEnergyInLC = 0.f;
0251           auto maxCPId = -1;
0252           for (auto& h : hit_find_in_CP->second) {
0253             auto shared_fraction = std::min(rhFraction, h.fraction);
0254             //We are in the case where there are caloParticles with simHits connected via detid with the recHit under study
0255             //So, from all layers clusters, find the recHits that are connected with a caloParticle and save/calculate the
0256             //energy of that caloParticle as the sum over all recHits of the recHits energy weighted
0257             //by the caloParticle's fraction related to that recHit.
0258             CPEnergyInMCL[h.clusterId] += shared_fraction * hit->energy();
0259             //Same but for layer clusters for the cell association per layer
0260             CPEnergyInLC[h.clusterId] += shared_fraction * hit->energy();
0261             //Here cPOnLayer[caloparticle][layer] described above is set
0262             //Here for multiClusters with matched recHit, the CP fraction times hit energy is added and saved
0263             cPOnLayer[h.clusterId][lcLayerId].multiClusterIdToEnergyAndScore[mcId].first +=
0264                 shared_fraction * hit->energy();
0265             cPOnLayer[h.clusterId][lcLayerId].multiClusterIdToEnergyAndScore[mcId].second = FLT_MAX;
0266             //cpsInMultiCluster[multicluster][CPids]
0267             //Connects a multiCluster with all related caloParticles
0268             cpsInMultiCluster[mcId].emplace_back(h.clusterId, FLT_MAX);
0269             //From all CaloParticles related to a layer cluster, we save id and energy of the caloParticle
0270             //that after simhit-rechit matching in layer has the maximum energy.
0271             if (shared_fraction > maxCPEnergyInLC) {
0272               //energy is used only here. cpid is saved for multiClusters
0273               maxCPEnergyInLC = CPEnergyInLC[h.clusterId];
0274               maxCPId = h.clusterId;
0275             }
0276           }
0277           //Keep in mind here maxCPId could be zero. So, below ask for negative not including zero to count noise.
0278           hitsToCaloParticleId[hitId] = maxCPId;
0279         }
0280 
0281       }  //end of loop through recHits of the layer cluster.
0282 
0283       //Loop through all recHits to count how many of them are noise and how many are matched.
0284       //In case of matched rechit-simhit, we count and save the number of recHits related to the maximum energy CaloParticle.
0285       for (auto c : hitsToCaloParticleId) {
0286         if (c < 0) {
0287           numberOfNoiseHitsInMCL++;
0288         } else {
0289           occurrencesCPinMCL[c]++;
0290         }
0291       }
0292 
0293       //Below from all maximum energy CaloParticles, we save the one with the largest amount
0294       //of related recHits.
0295       for (auto& c : occurrencesCPinMCL) {
0296         if (c.second > maxCPNumberOfHitsInMCL) {
0297           maxCPId_byNumberOfHits = c.first;
0298           maxCPNumberOfHitsInMCL = c.second;
0299         }
0300       }
0301 
0302       //Find the CaloParticle that has the maximum energy shared with the multiCluster under study.
0303       for (auto& c : CPEnergyInMCL) {
0304         if (c.second > maxEnergySharedMCLandCP) {
0305           maxCPId_byEnergy = c.first;
0306           maxEnergySharedMCLandCP = c.second;
0307         }
0308       }
0309       //The energy of the CaloParticle that found to have the maximum energy shared with the multiCluster under study.
0310       float totalCPEnergyFromLayerCP = 0.f;
0311       if (maxCPId_byEnergy >= 0) {
0312         //Loop through all layers
0313         for (unsigned int j = 0; j < layers_ * 2; ++j) {
0314           totalCPEnergyFromLayerCP = totalCPEnergyFromLayerCP + cPOnLayer[maxCPId_byEnergy][j].energy;
0315         }
0316         energyFractionOfCPinMCL = maxEnergySharedMCLandCP / totalCPEnergyFromLayerCP;
0317         if (mClusters[mcId].energy() > 0.f) {
0318           energyFractionOfMCLinCP = maxEnergySharedMCLandCP / mClusters[mcId].energy();
0319         }
0320       }
0321 
0322       LogDebug("MultiClusterAssociatorByEnergyScoreImpl") << std::setw(12) << "multiCluster"
0323                                                           << "\t" << std::setw(10) << "mulcl energy"
0324                                                           << "\t" << std::setw(5) << "nhits"
0325                                                           << "\t" << std::setw(12) << "noise hits"
0326                                                           << "\t" << std::setw(22) << "maxCPId_byNumberOfHits"
0327                                                           << "\t" << std::setw(8) << "nhitsCP"
0328                                                           << "\t" << std::setw(16) << "maxCPId_byEnergy"
0329                                                           << "\t" << std::setw(23) << "maxEnergySharedMCLandCP"
0330                                                           << "\t" << std::setw(22) << "totalCPEnergyFromAllLayerCP"
0331                                                           << "\t" << std::setw(22) << "energyFractionOfMCLinCP"
0332                                                           << "\t" << std::setw(25) << "energyFractionOfCPinMCL"
0333                                                           << "\t" << std::endl;
0334       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0335           << std::setw(12) << mcId << "\t" << std::setw(10) << mClusters[mcId].energy() << "\t" << std::setw(5)
0336           << numberOfHitsInMCL << "\t" << std::setw(12) << numberOfNoiseHitsInMCL << "\t" << std::setw(22)
0337           << maxCPId_byNumberOfHits << "\t" << std::setw(8) << maxCPNumberOfHitsInMCL << "\t" << std::setw(16)
0338           << maxCPId_byEnergy << "\t" << std::setw(23) << maxEnergySharedMCLandCP << "\t" << std::setw(22)
0339           << totalCPEnergyFromLayerCP << "\t" << std::setw(22) << energyFractionOfMCLinCP << "\t" << std::setw(25)
0340           << energyFractionOfCPinMCL << std::endl;
0341     }
0342   }  // end of loop through multiClusters
0343 
0344   // Update cpsInMultiCluster; compute the score MultiCluster-to-CaloParticle,
0345   // together with the returned AssociationMap
0346   for (unsigned int mcId = 0; mcId < nMultiClusters; ++mcId) {
0347     // find the unique caloParticles id contributing to the multilusters
0348     std::sort(cpsInMultiCluster[mcId].begin(), cpsInMultiCluster[mcId].end());
0349     auto last = std::unique(cpsInMultiCluster[mcId].begin(), cpsInMultiCluster[mcId].end());
0350     cpsInMultiCluster[mcId].erase(last, cpsInMultiCluster[mcId].end());
0351 
0352     const auto& hits_and_fractions = mClusters[mcId].hitsAndFractions();
0353     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0354     if (numberOfHitsInLC > 0) {
0355       if (mClusters[mcId].energy() == 0. && !cpsInMultiCluster[mcId].empty()) {
0356         //Loop through all CaloParticles contributing to multiCluster mcId.
0357         for (auto& cpPair : cpsInMultiCluster[mcId]) {
0358           //In case of a multiCluster with zero energy but related CaloParticles the score is set to 1.
0359           cpPair.second = 1.;
0360           LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0361               << "multiClusterId : \t " << mcId << "\t CP id : \t" << cpPair.first << "\t score \t " << cpPair.second
0362               << "\n";
0363         }
0364         continue;
0365       }
0366 
0367       // Compute the correct normalization
0368       float invMultiClusterEnergyWeight = 0.f;
0369       for (auto const& haf : mClusters[mcId].hitsAndFractions()) {
0370         invMultiClusterEnergyWeight +=
0371             (haf.second * hitMap_->at(haf.first)->energy()) * (haf.second * hitMap_->at(haf.first)->energy());
0372       }
0373       invMultiClusterEnergyWeight = 1.f / invMultiClusterEnergyWeight;
0374 
0375       for (unsigned int i = 0; i < numberOfHitsInLC; ++i) {
0376         DetId rh_detid = hits_and_fractions[i].first;
0377         float rhFraction = hits_and_fractions[i].second;
0378 
0379         bool hitWithNoCP = (detIdToCaloParticleId_Map.find(rh_detid) == detIdToCaloParticleId_Map.end());
0380 
0381         auto itcheck = hitMap_->find(rh_detid);
0382         const HGCRecHit* hit = itcheck->second;
0383         float hitEnergyWeight = hit->energy() * hit->energy();
0384 
0385         for (auto& cpPair : cpsInMultiCluster[mcId]) {
0386           unsigned int multiClusterId = cpPair.first;
0387           float cpFraction = 0.f;
0388           if (!hitWithNoCP) {
0389             auto findHitIt = std::find(detIdToCaloParticleId_Map[rh_detid].begin(),
0390                                        detIdToCaloParticleId_Map[rh_detid].end(),
0391                                        hgcal::detIdInfoInCluster{multiClusterId, 0.f});
0392             if (findHitIt != detIdToCaloParticleId_Map[rh_detid].end())
0393               cpFraction = findHitIt->fraction;
0394           }
0395           if (cpPair.second == FLT_MAX) {
0396             cpPair.second = 0.f;
0397           }
0398           cpPair.second +=
0399               (rhFraction - cpFraction) * (rhFraction - cpFraction) * hitEnergyWeight * invMultiClusterEnergyWeight;
0400         }
0401       }  // End of loop over Hits within a MultiCluster
0402 #ifdef EDM_ML_DEBUG
0403       //In case of a multiCluster with some energy but none related CaloParticles print some info.
0404       if (cpsInMultiCluster[mcId].empty())
0405         LogDebug("MultiClusterAssociatorByEnergyScoreImpl") << "multiCluster Id: \t" << mcId << "\tCP id:\t-1 "
0406                                                             << "\t score \t-1"
0407                                                             << "\n";
0408 #endif
0409     }
0410   }  // End of loop over MultiClusters
0411 
0412   // Compute the CaloParticle-To-MultiCluster score
0413   for (const auto& cpId : cPSelectedIndices) {
0414     for (unsigned int layerId = 0; layerId < layers_ * 2; ++layerId) {
0415       unsigned int CPNumberOfHits = cPOnLayer[cpId][layerId].hits_and_fractions.size();
0416       if (CPNumberOfHits == 0)
0417         continue;
0418 #ifdef EDM_ML_DEBUG
0419       int mcWithMaxEnergyInCP = -1;
0420       float maxEnergyMCLperlayerinCP = 0.f;
0421       float CPenergy = cPOnLayer[cpId][layerId].energy;
0422       float CPEnergyFractionInMCLperlayer = 0.f;
0423       for (auto& mc : cPOnLayer[cpId][layerId].multiClusterIdToEnergyAndScore) {
0424         if (mc.second.first > maxEnergyMCLperlayerinCP) {
0425           maxEnergyMCLperlayerinCP = mc.second.first;
0426           mcWithMaxEnergyInCP = mc.first;
0427         }
0428       }
0429       if (CPenergy > 0.f)
0430         CPEnergyFractionInMCLperlayer = maxEnergyMCLperlayerinCP / CPenergy;
0431 
0432       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0433           << std::setw(8) << "LayerId:\t" << std::setw(12) << "caloparticle\t" << std::setw(15) << "cp total energy\t"
0434           << std::setw(15) << "cpEnergyOnLayer\t" << std::setw(14) << "CPNhitsOnLayer\t" << std::setw(18)
0435           << "mcWithMaxEnergyInCP\t" << std::setw(15) << "maxEnergyMCLinCP\t" << std::setw(20)
0436           << "CPEnergyFractionInMCL"
0437           << "\n";
0438       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0439           << std::setw(8) << layerId << "\t" << std::setw(12) << cpId << "\t" << std::setw(15)
0440           << caloParticles[cpId].energy() << "\t" << std::setw(15) << CPenergy << "\t" << std::setw(14)
0441           << CPNumberOfHits << "\t" << std::setw(18) << mcWithMaxEnergyInCP << "\t" << std::setw(15)
0442           << maxEnergyMCLperlayerinCP << "\t" << std::setw(20) << CPEnergyFractionInMCLperlayer << "\n";
0443 #endif
0444 
0445       for (unsigned int i = 0; i < CPNumberOfHits; ++i) {
0446         auto& cp_hitDetId = cPOnLayer[cpId][layerId].hits_and_fractions[i].first;
0447         auto& cpFraction = cPOnLayer[cpId][layerId].hits_and_fractions[i].second;
0448 
0449         bool hitWithNoMCL = false;
0450         if (cpFraction == 0.f)
0451           continue;  //hopefully this should never happen
0452         auto hit_find_in_MCL = detIdToMultiClusterId_Map.find(cp_hitDetId);
0453         if (hit_find_in_MCL == detIdToMultiClusterId_Map.end())
0454           hitWithNoMCL = true;
0455         auto itcheck = hitMap_->find(cp_hitDetId);
0456         const HGCRecHit* hit = itcheck->second;
0457         float hitEnergyWeight = hit->energy() * hit->energy();
0458         for (auto& mcPair : cPOnLayer[cpId][layerId].multiClusterIdToEnergyAndScore) {
0459           unsigned int multiClusterId = mcPair.first;
0460           float mcFraction = 0.f;
0461 
0462           if (!hitWithNoMCL) {
0463             auto findHitIt = std::find(detIdToMultiClusterId_Map[cp_hitDetId].begin(),
0464                                        detIdToMultiClusterId_Map[cp_hitDetId].end(),
0465                                        hgcal::detIdInfoInMultiCluster{multiClusterId, 0, 0.f});
0466             if (findHitIt != detIdToMultiClusterId_Map[cp_hitDetId].end())
0467               mcFraction = findHitIt->fraction;
0468           }
0469           //Observe here that we do not divide as before by the layer cluster energy weight. We should sum first
0470           //over all layers and divide with the total CP energy over all layers.
0471           if (mcPair.second.second == FLT_MAX) {
0472             mcPair.second.second = 0.f;
0473           }
0474           mcPair.second.second += (mcFraction - cpFraction) * (mcFraction - cpFraction) * hitEnergyWeight;
0475 #ifdef EDM_ML_DEBUG
0476           LogDebug("HGCalValidator") << "multiClusterId:\t" << multiClusterId << "\tmcfraction,cpfraction:\t"
0477                                      << mcFraction << ", " << cpFraction << "\thitEnergyWeight:\t" << hitEnergyWeight
0478                                      << "\tcurrent score numerator:\t" << mcPair.second.second << "\n";
0479 #endif
0480         }  // End of loop over MultiClusters linked to hits of this CaloParticle
0481       }    // End of loop over hits of CaloParticle on a Layer
0482 #ifdef EDM_ML_DEBUG
0483       if (cPOnLayer[cpId][layerId].multiClusterIdToEnergyAndScore.empty())
0484         LogDebug("HGCalValidator") << "CP Id: \t" << cpId << "\t MCL id:\t-1 "
0485                                    << "\t layer \t " << layerId << " Sub score in \t -1"
0486                                    << "\n";
0487 
0488 #endif
0489     }
0490   }
0491   return {cpsInMultiCluster, cPOnLayer};
0492 }
0493 
0494 hgcal::RecoToSimCollectionWithMultiClusters MultiClusterAssociatorByEnergyScoreImpl::associateRecoToSim(
0495     const edm::Handle<reco::HGCalMultiClusterCollection>& mCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0496   hgcal::RecoToSimCollectionWithMultiClusters returnValue(productGetter_);
0497   const auto& links = makeConnections(mCCH, cPCH);
0498 
0499   const auto& cpsInMultiCluster = std::get<0>(links);
0500   for (size_t mcId = 0; mcId < cpsInMultiCluster.size(); ++mcId) {
0501     for (auto& cpPair : cpsInMultiCluster[mcId]) {
0502       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0503           << "multiCluster Id: \t" << mcId << "\t CP id: \t" << cpPair.first << "\t score \t" << cpPair.second << "\n";
0504       // Fill AssociationMap
0505       returnValue.insert(edm::Ref<reco::HGCalMultiClusterCollection>(mCCH, mcId),  // Ref to MC
0506                          std::make_pair(edm::Ref<CaloParticleCollection>(cPCH, cpPair.first),
0507                                         cpPair.second)  // Pair <Ref to CP, score>
0508       );
0509     }
0510   }
0511   return returnValue;
0512 }
0513 
0514 hgcal::SimToRecoCollectionWithMultiClusters MultiClusterAssociatorByEnergyScoreImpl::associateSimToReco(
0515     const edm::Handle<reco::HGCalMultiClusterCollection>& mCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0516   hgcal::SimToRecoCollectionWithMultiClusters returnValue(productGetter_);
0517   const auto& links = makeConnections(mCCH, cPCH);
0518   const auto& cPOnLayer = std::get<1>(links);
0519   for (size_t cpId = 0; cpId < cPOnLayer.size(); ++cpId) {
0520     for (size_t layerId = 0; layerId < cPOnLayer[cpId].size(); ++layerId) {
0521       for (auto& mcPair : cPOnLayer[cpId][layerId].multiClusterIdToEnergyAndScore) {
0522         returnValue.insert(
0523             edm::Ref<CaloParticleCollection>(cPCH, cpId),                                    // Ref to CP
0524             std::make_pair(edm::Ref<reco::HGCalMultiClusterCollection>(mCCH, mcPair.first),  // Pair <Ref to MC,
0525                            std::make_pair(mcPair.second.first, mcPair.second.second))        // pair <energy, score> >
0526         );
0527       }
0528     }
0529   }
0530   return returnValue;
0531 }