HGCalAssociatorProducers/plugins/TSToSCAssociatorByEnergyScoreImpl.cc

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