HGCalAssociatorProducers/plugins/MultiClusterAssociatorByEnergyScoreImpl.cc

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 numberOfHitsInMCL = 0;
0184
0185       //number of hits related to that cluster
0186       unsigned int numberOfHitsInLC = hits_and_fractions.size();
0187       numberOfHitsInMCL += numberOfHitsInLC;
0188       std::unordered_map<unsigned, float> CPEnergyInLC;
0189
0190       //hitsToCaloParticleId is a vector of ints, one for each recHit of the
0191       //layer cluster under study. If negative, there is no simHit from any CaloParticle related.
0192       //If positive, at least one CaloParticle has been found with matched simHit.
0193       //In more detail:
0194       // 1. hitsToCaloParticleId[hitId] = -3
0195       //    TN:  These represent Halo Cells(N) that have not been
0196       //    assigned to any CaloParticle (hence the T).
0197       // 2. hitsToCaloParticleId[hitId] = -2
0198       //    FN: There represent Halo Cells(N) that have been assigned
0199       //    to a CaloParticle (hence the F, since those should have not been marked as halo)
0200       // 3. hitsToCaloParticleId[hitId] = -1
0201       //    FP: These represent Real Cells(P) that have not been
0202       //    assigned to any CaloParticle (hence the F, since these are fakes)
0203       // 4. hitsToCaloParticleId[hitId] >= 0
0204       //    TP There represent Real Cells(P) that have been assigned
0205       //    to a CaloParticle (hence the T)
0206
0207       std::vector<int> hitsToCaloParticleId(numberOfHitsInLC);
0208       //det id of the first hit just to make the lcLayerId variable
0209       //which maps the layers in -z: 0->51 and in +z: 52->103
0210       const auto firstHitDetId = hits_and_fractions[0].first;
0211       int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) +
0212                       layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
0213
0214       //Loop though the hits of the layer cluster under study
0215       for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0216         const auto rh_detid = hits_and_fractions[hitId].first;
0217         const auto rhFraction = hits_and_fractions[hitId].second;
0218
0219         //Since the hit is belonging to the layer cluster, it must also be in the recHits map.
0220         const auto itcheck = hitMap_->find(rh_detid);
0221         const auto hit = itcheck->second;
0222
0223         //Make a map that will connect a detid (that belongs to a recHit of the layer cluster under study,
0224         //no need to save others) with:
0225         //1. the layer clusters that have recHits in that detid
0226         //2. the fraction of the recHit of each layer cluster that contributes to that detid.
0227         //So, something like:
0228         //detid: (layer cluster 1, hit fraction) , (layer cluster 2, hit fraction), (layer cluster 3, hit fraction) ...
0229         //here comparing with the caloParticle map above
0230         auto hit_find_in_LC = detIdToMultiClusterId_Map.find(rh_detid);
0231         if (hit_find_in_LC == detIdToMultiClusterId_Map.end()) {
0232           detIdToMultiClusterId_Map[rh_detid] = std::vector<hgcal::detIdInfoInMultiCluster>();
0233         }
0234         detIdToMultiClusterId_Map[rh_detid].emplace_back(hgcal::detIdInfoInMultiCluster{mcId, mcId, rhFraction});
0235
0236         // Check whether the recHit of the layer cluster under study has a sim hit in the same cell
0237         auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
0238
0239         // If the fraction is zero or the hit does not belong to any calo
0240         // particle, set the caloParticleId for the hit to -1 and this will
0241         // contribute to the number of noise hits
0242         if (rhFraction == 0.) {  // this could be a real hit that has been marked as halo
0243           hitsToCaloParticleId[hitId] = -2;
0244         }
0245         if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) {
0246           hitsToCaloParticleId[hitId] -= 1;
0247         } else {
0248           auto maxCPEnergyInLC = 0.f;
0249           auto maxCPId = -1;
0250           for (auto& h : hit_find_in_CP->second) {
0251             auto shared_fraction = std::min(rhFraction, h.fraction);
0252             //We are in the case where there are caloParticles with simHits connected via detid with the recHit under study
0253             //So, from all layers clusters, find the recHits that are connected with a caloParticle and save/calculate the
0254             //energy of that caloParticle as the sum over all recHits of the recHits energy weighted
0255             //by the caloParticle's fraction related to that recHit.
0256             CPEnergyInMCL[h.clusterId] += shared_fraction * hit->energy();
0257             //Same but for layer clusters for the cell association per layer
0258             CPEnergyInLC[h.clusterId] += shared_fraction * hit->energy();
0259             //Here cPOnLayer[caloparticle][layer] described above is set
0260             //Here for multiClusters with matched recHit, the CP fraction times hit energy is added and saved
0261             cPOnLayer[h.clusterId][lcLayerId].multiClusterIdToEnergyAndScore[mcId].first +=
0262                 shared_fraction * hit->energy();
0263             cPOnLayer[h.clusterId][lcLayerId].multiClusterIdToEnergyAndScore[mcId].second = FLT_MAX;
0264             //cpsInMultiCluster[multicluster][CPids]
0265             //Connects a multiCluster with all related caloParticles
0266             cpsInMultiCluster[mcId].emplace_back(h.clusterId, FLT_MAX);
0267             //From all CaloParticles related to a layer cluster, we save id and energy of the caloParticle
0268             //that after simhit-rechit matching in layer has the maximum energy.
0269             if (shared_fraction > maxCPEnergyInLC) {
0270               //energy is used only here. cpid is saved for multiClusters
0271               maxCPEnergyInLC = CPEnergyInLC[h.clusterId];
0272               maxCPId = h.clusterId;
0273             }
0274           }
0275           //Keep in mind here maxCPId could be zero. So, below ask for negative not including zero to count noise.
0276           hitsToCaloParticleId[hitId] = maxCPId;
0277         }
0278
0279       }  //end of loop through recHits of the layer cluster.
0280
0281       //Loop through all recHits to count how many of them are noise and how many are matched.
0282       //In case of matched rechit-simhit, we count and save the number of recHits related to the maximum energy CaloParticle.
0283       for (auto c : hitsToCaloParticleId) {
0284         if (c < 0) {
0285           numberOfNoiseHitsInMCL++;
0286         } else {
0287           occurrencesCPinMCL[c]++;
0288         }
0289       }
0290
0291       //Below from all maximum energy CaloParticles, we save the one with the largest amount
0292       //of related recHits.
0293       for (auto& c : occurrencesCPinMCL) {
0294         if (c.second > maxCPNumberOfHitsInMCL) {
0295           maxCPId_byNumberOfHits = c.first;
0296           maxCPNumberOfHitsInMCL = c.second;
0297         }
0298       }
0299
0300       //Find the CaloParticle that has the maximum energy shared with the multiCluster under study.
0301       for (auto& c : CPEnergyInMCL) {
0302         if (c.second > maxEnergySharedMCLandCP) {
0303           maxCPId_byEnergy = c.first;
0304           maxEnergySharedMCLandCP = c.second;
0305         }
0306       }
0307       //The energy of the CaloParticle that found to have the maximum energy shared with the multiCluster under study.
0308       float totalCPEnergyFromLayerCP = 0.f;
0309       if (maxCPId_byEnergy >= 0) {
0310         //Loop through all layers
0311         for (unsigned int j = 0; j < layers_ * 2; ++j) {
0312           totalCPEnergyFromLayerCP = totalCPEnergyFromLayerCP + cPOnLayer[maxCPId_byEnergy][j].energy;
0313         }
0314         energyFractionOfCPinMCL = maxEnergySharedMCLandCP / totalCPEnergyFromLayerCP;
0315         if (mClusters[mcId].energy() > 0.f) {
0316           energyFractionOfMCLinCP = maxEnergySharedMCLandCP / mClusters[mcId].energy();
0317         }
0318       }
0319
0320       LogDebug("MultiClusterAssociatorByEnergyScoreImpl") << std::setw(12) << "multiCluster"
0321                                                           << "\t" << std::setw(10) << "mulcl energy"
0322                                                           << "\t" << std::setw(5) << "nhits"
0323                                                           << "\t" << std::setw(12) << "noise hits"
0324                                                           << "\t" << std::setw(22) << "maxCPId_byNumberOfHits"
0325                                                           << "\t" << std::setw(8) << "nhitsCP"
0326                                                           << "\t" << std::setw(16) << "maxCPId_byEnergy"
0327                                                           << "\t" << std::setw(23) << "maxEnergySharedMCLandCP"
0328                                                           << "\t" << std::setw(22) << "totalCPEnergyFromAllLayerCP"
0329                                                           << "\t" << std::setw(22) << "energyFractionOfMCLinCP"
0330                                                           << "\t" << std::setw(25) << "energyFractionOfCPinMCL"
0331                                                           << "\t" << std::endl;
0332       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0333           << std::setw(12) << mcId << "\t" << std::setw(10) << mClusters[mcId].energy() << "\t" << std::setw(5)
0334           << numberOfHitsInMCL << "\t" << std::setw(12) << numberOfNoiseHitsInMCL << "\t" << std::setw(22)
0335           << maxCPId_byNumberOfHits << "\t" << std::setw(8) << maxCPNumberOfHitsInMCL << "\t" << std::setw(16)
0336           << maxCPId_byEnergy << "\t" << std::setw(23) << maxEnergySharedMCLandCP << "\t" << std::setw(22)
0337           << totalCPEnergyFromLayerCP << "\t" << std::setw(22) << energyFractionOfMCLinCP << "\t" << std::setw(25)
0338           << energyFractionOfCPinMCL << std::endl;
0339     }
0340   }  // end of loop through multiClusters
0341
0342   // Update cpsInMultiCluster; compute the score MultiCluster-to-CaloParticle,
0343   // together with the returned AssociationMap
0344   for (unsigned int mcId = 0; mcId < nMultiClusters; ++mcId) {
0345     // find the unique caloParticles id contributing to the multilusters
0346     std::sort(cpsInMultiCluster[mcId].begin(), cpsInMultiCluster[mcId].end());
0347     auto last = std::unique(cpsInMultiCluster[mcId].begin(), cpsInMultiCluster[mcId].end());
0348     cpsInMultiCluster[mcId].erase(last, cpsInMultiCluster[mcId].end());
0349
0350     const auto& hits_and_fractions = mClusters[mcId].hitsAndFractions();
0351     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0352     if (numberOfHitsInLC > 0) {
0353       if (mClusters[mcId].energy() == 0. && !cpsInMultiCluster[mcId].empty()) {
0354         //Loop through all CaloParticles contributing to multiCluster mcId.
0355         for (auto& cpPair : cpsInMultiCluster[mcId]) {
0356           //In case of a multiCluster with zero energy but related CaloParticles the score is set to 1.
0357           cpPair.second = 1.;
0358           LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0359               << "multiClusterId : \t " << mcId << "\t CP id : \t" << cpPair.first << "\t score \t " << cpPair.second
0360               << "\n";
0361         }
0362         continue;
0363       }
0364
0365       // Compute the correct normalization
0366       float invMultiClusterEnergyWeight = 0.f;
0367       for (auto const& haf : mClusters[mcId].hitsAndFractions()) {
0368         invMultiClusterEnergyWeight +=
0369             (haf.second * hitMap_->at(haf.first)->energy()) * (haf.second * hitMap_->at(haf.first)->energy());
0370       }
0371       invMultiClusterEnergyWeight = 1.f / invMultiClusterEnergyWeight;
0372
0373       for (unsigned int i = 0; i < numberOfHitsInLC; ++i) {
0374         DetId rh_detid = hits_and_fractions[i].first;
0375         float rhFraction = hits_and_fractions[i].second;
0376
0377         bool hitWithNoCP = (detIdToCaloParticleId_Map.find(rh_detid) == detIdToCaloParticleId_Map.end());
0378
0379         auto itcheck = hitMap_->find(rh_detid);
0380         const HGCRecHit* hit = itcheck->second;
0381         float hitEnergyWeight = hit->energy() * hit->energy();
0382
0383         for (auto& cpPair : cpsInMultiCluster[mcId]) {
0384           unsigned int multiClusterId = cpPair.first;
0385           float cpFraction = 0.f;
0386           if (!hitWithNoCP) {
0387             auto findHitIt = std::find(detIdToCaloParticleId_Map[rh_detid].begin(),
0388                                        detIdToCaloParticleId_Map[rh_detid].end(),
0389                                        hgcal::detIdInfoInCluster{multiClusterId, 0.f});
0390             if (findHitIt != detIdToCaloParticleId_Map[rh_detid].end())
0391               cpFraction = findHitIt->fraction;
0392           }
0393           if (cpPair.second == FLT_MAX) {
0394             cpPair.second = 0.f;
0395           }
0396           cpPair.second +=
0397               (rhFraction - cpFraction) * (rhFraction - cpFraction) * hitEnergyWeight * invMultiClusterEnergyWeight;
0398         }
0399       }  // End of loop over Hits within a MultiCluster
0400 #ifdef EDM_ML_DEBUG
0401       //In case of a multiCluster with some energy but none related CaloParticles print some info.
0402       if (cpsInMultiCluster[mcId].empty())
0403         LogDebug("MultiClusterAssociatorByEnergyScoreImpl") << "multiCluster Id: \t" << mcId << "\tCP id:\t-1 "
0404                                                             << "\t score \t-1"
0405                                                             << "\n";
0406 #endif
0407     }
0408   }  // End of loop over MultiClusters
0409
0410   // Compute the CaloParticle-To-MultiCluster score
0411   for (const auto& cpId : cPSelectedIndices) {
0412     for (unsigned int layerId = 0; layerId < layers_ * 2; ++layerId) {
0413       unsigned int CPNumberOfHits = cPOnLayer[cpId][layerId].hits_and_fractions.size();
0414       if (CPNumberOfHits == 0)
0415         continue;
0416 #ifdef EDM_ML_DEBUG
0417       int mcWithMaxEnergyInCP = -1;
0418       float maxEnergyMCLperlayerinCP = 0.f;
0419       float CPenergy = cPOnLayer[cpId][layerId].energy;
0420       float CPEnergyFractionInMCLperlayer = 0.f;
0421       for (auto& mc : cPOnLayer[cpId][layerId].multiClusterIdToEnergyAndScore) {
0422         if (mc.second.first > maxEnergyMCLperlayerinCP) {
0423           maxEnergyMCLperlayerinCP = mc.second.first;
0424           mcWithMaxEnergyInCP = mc.first;
0425         }
0426       }
0427       if (CPenergy > 0.f)
0428         CPEnergyFractionInMCLperlayer = maxEnergyMCLperlayerinCP / CPenergy;
0429
0430       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0431           << std::setw(8) << "LayerId:\t" << std::setw(12) << "caloparticle\t" << std::setw(15) << "cp total energy\t"
0432           << std::setw(15) << "cpEnergyOnLayer\t" << std::setw(14) << "CPNhitsOnLayer\t" << std::setw(18)
0433           << "mcWithMaxEnergyInCP\t" << std::setw(15) << "maxEnergyMCLinCP\t" << std::setw(20)
0434           << "CPEnergyFractionInMCL"
0435           << "\n";
0436       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0437           << std::setw(8) << layerId << "\t" << std::setw(12) << cpId << "\t" << std::setw(15)
0438           << caloParticles[cpId].energy() << "\t" << std::setw(15) << CPenergy << "\t" << std::setw(14)
0439           << CPNumberOfHits << "\t" << std::setw(18) << mcWithMaxEnergyInCP << "\t" << std::setw(15)
0440           << maxEnergyMCLperlayerinCP << "\t" << std::setw(20) << CPEnergyFractionInMCLperlayer << "\n";
0441 #endif
0442
0443       for (unsigned int i = 0; i < CPNumberOfHits; ++i) {
0444         auto& cp_hitDetId = cPOnLayer[cpId][layerId].hits_and_fractions[i].first;
0445         auto& cpFraction = cPOnLayer[cpId][layerId].hits_and_fractions[i].second;
0446
0447         bool hitWithNoMCL = false;
0448         if (cpFraction == 0.f)
0449           continue;  //hopefully this should never happen
0450         auto hit_find_in_MCL = detIdToMultiClusterId_Map.find(cp_hitDetId);
0451         if (hit_find_in_MCL == detIdToMultiClusterId_Map.end())
0452           hitWithNoMCL = true;
0453         auto itcheck = hitMap_->find(cp_hitDetId);
0454         const HGCRecHit* hit = itcheck->second;
0455         float hitEnergyWeight = hit->energy() * hit->energy();
0456         for (auto& mcPair : cPOnLayer[cpId][layerId].multiClusterIdToEnergyAndScore) {
0457           unsigned int multiClusterId = mcPair.first;
0458           float mcFraction = 0.f;
0459
0460           if (!hitWithNoMCL) {
0461             auto findHitIt = std::find(detIdToMultiClusterId_Map[cp_hitDetId].begin(),
0462                                        detIdToMultiClusterId_Map[cp_hitDetId].end(),
0463                                        hgcal::detIdInfoInMultiCluster{multiClusterId, 0, 0.f});
0464             if (findHitIt != detIdToMultiClusterId_Map[cp_hitDetId].end())
0465               mcFraction = findHitIt->fraction;
0466           }
0467           //Observe here that we do not divide as before by the layer cluster energy weight. We should sum first
0468           //over all layers and divide with the total CP energy over all layers.
0469           if (mcPair.second.second == FLT_MAX) {
0470             mcPair.second.second = 0.f;
0471           }
0472           mcPair.second.second += (mcFraction - cpFraction) * (mcFraction - cpFraction) * hitEnergyWeight;
0473 #ifdef EDM_ML_DEBUG
0474           LogDebug("HGCalValidator") << "multiClusterId:\t" << multiClusterId << "\tmcfraction,cpfraction:\t"
0475                                      << mcFraction << ", " << cpFraction << "\thitEnergyWeight:\t" << hitEnergyWeight
0476                                      << "\tcurrent score numerator:\t" << mcPair.second.second << "\n";
0477 #endif
0478         }  // End of loop over MultiClusters linked to hits of this CaloParticle
0479       }    // End of loop over hits of CaloParticle on a Layer
0480 #ifdef EDM_ML_DEBUG
0481       if (cPOnLayer[cpId][layerId].multiClusterIdToEnergyAndScore.empty())
0482         LogDebug("HGCalValidator") << "CP Id: \t" << cpId << "\t MCL id:\t-1 "
0483                                    << "\t layer \t " << layerId << " Sub score in \t -1"
0484                                    << "\n";
0485
0486 #endif
0487     }
0488   }
0489   return {cpsInMultiCluster, cPOnLayer};
0490 }
0491
0492 hgcal::RecoToSimCollectionWithMultiClusters MultiClusterAssociatorByEnergyScoreImpl::associateRecoToSim(
0493     const edm::Handle<reco::HGCalMultiClusterCollection>& mCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0494   hgcal::RecoToSimCollectionWithMultiClusters returnValue(productGetter_);
0495   const auto& links = makeConnections(mCCH, cPCH);
0496
0497   const auto& cpsInMultiCluster = std::get<0>(links);
0498   for (size_t mcId = 0; mcId < cpsInMultiCluster.size(); ++mcId) {
0499     for (auto& cpPair : cpsInMultiCluster[mcId]) {
0500       LogDebug("MultiClusterAssociatorByEnergyScoreImpl")
0501           << "multiCluster Id: \t" << mcId << "\t CP id: \t" << cpPair.first << "\t score \t" << cpPair.second << "\n";
0502       // Fill AssociationMap
0503       returnValue.insert(edm::Ref<reco::HGCalMultiClusterCollection>(mCCH, mcId),  // Ref to MC
0504                          std::make_pair(edm::Ref<CaloParticleCollection>(cPCH, cpPair.first),
0505                                         cpPair.second)  // Pair <Ref to CP, score>
0506       );
0507     }
0508   }
0509   return returnValue;
0510 }
0511
0512 hgcal::SimToRecoCollectionWithMultiClusters MultiClusterAssociatorByEnergyScoreImpl::associateSimToReco(
0513     const edm::Handle<reco::HGCalMultiClusterCollection>& mCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0514   hgcal::SimToRecoCollectionWithMultiClusters returnValue(productGetter_);
0515   const auto& links = makeConnections(mCCH, cPCH);
0516   const auto& cPOnLayer = std::get<1>(links);
0517   for (size_t cpId = 0; cpId < cPOnLayer.size(); ++cpId) {
0518     for (size_t layerId = 0; layerId < cPOnLayer[cpId].size(); ++layerId) {
0519       for (auto& mcPair : cPOnLayer[cpId][layerId].multiClusterIdToEnergyAndScore) {
0520         returnValue.insert(
0521             edm::Ref<CaloParticleCollection>(cPCH, cpId),                                    // Ref to CP
0522             std::make_pair(edm::Ref<reco::HGCalMultiClusterCollection>(mCCH, mcPair.first),  // Pair <Ref to MC,
0523                            std::make_pair(mcPair.second.first, mcPair.second.second))        // pair <energy, score> >
0524         );
0525       }
0526     }
0527   }
0528   return returnValue;
0529 }