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