HGCalAssociatorProducers/plugins/TSToSimTSAssociatorByEnergyScoreImpl.cc

0001 #include "TSToSimTSAssociatorByEnergyScoreImpl.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 TSToSimTSAssociatorByEnergyScoreImpl::TSToSimTSAssociatorByEnergyScoreImpl(
0006     edm::EDProductGetter const& productGetter,
0007     bool hardScatterOnly,
0008     std::shared_ptr<hgcal::RecHitTools> recHitTools,
0009     const std::unordered_map<DetId, const unsigned int>* hitMap,
0010     std::vector<const HGCRecHit*>& hits)
0011     : hardScatterOnly_(hardScatterOnly),
0012       recHitTools_(recHitTools),
0013       hitMap_(hitMap),
0014       hits_(hits),
0015       productGetter_(&productGetter) {
0016   layers_ = recHitTools_->lastLayerBH();
0017 }
0018
0019 ticl::association TSToSimTSAssociatorByEnergyScoreImpl::makeConnections(
0020     const edm::Handle<ticl::TracksterCollection>& tCH,
0021     const edm::Handle<reco::CaloClusterCollection>& lCCH,
0022     const edm::Handle<ticl::TracksterCollection>& sTCH) const {
0023   // Get collections
0024   const auto& tracksters = *tCH.product();
0025   const auto& layerClusters = *lCCH.product();
0026   const auto& simTracksters = *sTCH.product();
0027   auto nTracksters = tracksters.size();
0028
0029   //There shouldn't be any SimTrackster without vertices, but maybe they will be added later.
0030   auto nSimTracksters = simTracksters.size();
0031   std::vector<size_t> sTIndices;
0032   for (unsigned int stId = 0; stId < nSimTracksters; ++stId) {
0033     if (simTracksters[stId].vertices().empty()) {
0034       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0035           << "Excluding SimTrackster " << stId << " witH no vertices!" << std::endl;
0036       continue;
0037     }
0038     sTIndices.emplace_back(stId);
0039   }
0040   nSimTracksters = sTIndices.size();
0041
0042   // Initialize tssInSimTrackster. It contains the simTracksterOnLayer structure for all CaloParticles in each layer and
0043   // among other the information to compute the SimTrackster-To-Trackster score. It is one of the two objects that
0044   // build the output of the makeConnections function.
0045   // tssInSimTrackster[stId]
0046   ticl::simTracksterToTrackster tssInSimTrackster;
0047   tssInSimTrackster.resize(nSimTracksters);
0048   for (unsigned int i = 0; i < nSimTracksters; ++i) {
0049     tssInSimTrackster[i].simTracksterId = i;
0050     tssInSimTrackster[i].energy = 0.f;
0051     tssInSimTrackster[i].lcs_and_fractions.clear();
0052   }
0053
0054   // Fill lcToSimTracksterId_Map and update tssInSimTrackster
0055   // The lcToSimTracksterId_Map is used to connect a LayerCluster, via its id (key), with all the SimTracksters that
0056   // contributed to that LayerCluster by storing the SimTrackster id and the fraction of the LayerCluster's energy
0057   // in which the SimTrackster contributed.
0058   std::unordered_map<int, std::vector<ticl::lcInfoInTrackster>> lcToSimTracksterId_Map;
0059   for (const auto& stId : sTIndices) {
0060     const auto& lcs = simTracksters[stId].vertices();
0061     int lcInSimTst = 0;
0062     for (const auto& lcId : lcs) {
0063       const auto fraction = 1.f / simTracksters[stId].vertex_multiplicity(lcInSimTst++);
0064
0065       const auto lc_find_it = lcToSimTracksterId_Map.find(lcId);
0066       if (lc_find_it == lcToSimTracksterId_Map.end()) {
0067         lcToSimTracksterId_Map[lcId] = std::vector<ticl::lcInfoInTrackster>();
0068       }
0069       lcToSimTracksterId_Map[lcId].emplace_back(stId, fraction);
0070
0071       tssInSimTrackster[stId].energy += fraction * layerClusters[lcId].energy();
0072       tssInSimTrackster[stId].lcs_and_fractions.emplace_back(lcId, fraction);
0073     }
0074   }  // end of loop over SimTracksters
0075
0076 #ifdef EDM_ML_DEBUG
0077   LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0078       << "tssInSimTrackster INFO (Only SimTrackster filled at the moment)" << std::endl;
0079   for (size_t st = 0; st < tssInSimTrackster.size(); ++st) {
0080     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "For SimTrackster Idx: " << st << " we have: " << std::endl;
0081     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0082         << "\tSimTracksterIdx:\t" << tssInSimTrackster[st].simTracksterId << std::endl;
0083     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "\tEnergy:\t" << tssInSimTrackster[st].energy << std::endl;
0084     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "\t# of clusters:\t" << layerClusters.size() << std::endl;
0085     double tot_energy = 0.;
0086     for (auto const& lcaf : tssInSimTrackster[st].lcs_and_fractions) {
0087       const auto lcEn = lcaf.second * layerClusters[lcaf.first].energy();
0088       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0089           << "\tLC/fraction/energy: " << (uint32_t)lcaf.first << "/" << lcaf.second << "/" << lcEn << std::endl;
0090       tot_energy += lcEn;
0091     }
0092     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "\tTot Sum lcaf: " << tot_energy << std::endl;
0093     for (auto const& ts : tssInSimTrackster[st].tracksterIdToEnergyAndScore) {
0094       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0095           << "\ttsIdx/energy/score: " << ts.first << "/" << ts.second.first << "/" << ts.second.second << std::endl;
0096     }
0097   }
0098
0099   LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "lcToSimTracksterId_Map INFO" << std::endl;
0100   for (auto const& lc : lcToSimTracksterId_Map) {
0101     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0102         << "For lcId: " << (uint32_t)lc.first
0103         << " we have found the following connections with SimTracksters:" << std::endl;
0104     for (auto const& st : lc.second) {
0105       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0106           << "\tSimTrackster Id: " << st.clusterId << " with fraction: " << st.fraction
0107           << " and energy: " << st.fraction * layerClusters[lc.first].energy() << std::endl;
0108     }
0109   }
0110 #endif
0111
0112   // Fill stsInTrackster and update tssInSimTrackster
0113   // this contains the ids of the simTracksters contributing with at least one
0114   // hit to the Trackster. To be returned since this contains the information
0115   // to compute the Trackster-To-SimTrackster score.
0116   ticl::tracksterToSimTrackster stsInTrackster;  // tsId->(stId,score)
0117   stsInTrackster.resize(nTracksters);
0118
0119   for (unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
0120     for (unsigned int i = 0; i < tracksters[tsId].vertices().size(); ++i) {
0121       const auto lcId = tracksters[tsId].vertices(i);
0122       const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(i);
0123
0124       const auto lc_find_in_ST = lcToSimTracksterId_Map.find(lcId);
0125
0126       if (lc_find_in_ST != lcToSimTracksterId_Map.end()) {
0127         //Loops through all the simTracksters that contain the layer cluster under study
0128         //Here is time to update the tssInSimTrackster and connect the SimTrackster with all
0129         //the Tracksters that have the current layer cluster matched.
0130         for (const auto& st : lc_find_in_ST->second) {
0131           //tssInSimTrackster[simTracksterId][layerclusterId]-> (energy,score)
0132           //ST_i - > TS_j, TS_k, ...
0133           tssInSimTrackster[st.clusterId].tracksterIdToEnergyAndScore[tsId].first +=
0134               lcFractionInTs * st.fraction * layerClusters[lcId].energy();
0135           //TS_i -> ST_j, ST_k, ...
0136           stsInTrackster[tsId].emplace_back(st.clusterId, std::make_pair(0.f, 0.f));
0137         }
0138       }
0139     }  // End loop over LayerClusters in Trackster
0140   }  // End of loop over Tracksters
0141
0142 #ifdef EDM_ML_DEBUG
0143   for (unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
0144     for (const auto& lcId : tracksters[tsId].vertices()) {
0145       const auto& hits_and_fractions = layerClusters[lcId].hitsAndFractions();
0146       unsigned int numberOfHitsInLC = hits_and_fractions.size();
0147
0148       // This vector will store, for each hit in the Layercluster, the index of
0149       // the SimTrackster that contributed the most, in terms of energy, to it.
0150       // Special values are:
0151       //
0152       // -2  --> the reconstruction fraction of the RecHit is 0 (used in the past to monitor Halo Hits)
0153       // -3  --> same as before with the added condition that no SimTrackster has been linked to this RecHit
0154       // -1  --> the reco fraction is >0, but no SimTrackster has been linked to it
0155       // >=0 --> index of the linked SimTrackster
0156       std::vector<int> hitsToSimTracksterId(numberOfHitsInLC);
0157       // This will store the index of the SimTrackster linked to the LayerCluster that has the largest number of hits in common.
0158       int maxSTId_byNumberOfHits = -1;
0159       // This will store the maximum number of shared hits between a LayerCluster and a SimTrackster
0160       unsigned int maxSTNumberOfHitsInLC = 0;
0161       // This will store the index of the SimTrackster linked to the LayerCluster that has the largest energy in common.
0162       int maxSTId_byEnergy = -1;
0163       // This will store the maximum number of shared energy between a LayerCluster and a SimTrackster
0164       float maxEnergySharedLCandST = 0.f;
0165       // This will store the fraction of the LayerCluster energy shared with the best(energy) SimTrackster: e_shared/lc_energy
0166       float energyFractionOfLCinST = 0.f;
0167       // This will store the fraction of the SimTrackster energy shared with the Trackster: e_shared/sc_energy
0168       float energyFractionOfSTinLC = 0.f;
0169       std::unordered_map<unsigned, unsigned> occurrencesSTinLC;
0170       unsigned int numberOfNoiseHitsInLC = 0;
0171       std::unordered_map<unsigned, float> STEnergyInLC;
0172
0173       const auto lc_find_in_ST = lcToSimTracksterId_Map.find(lcId);
0174       for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0175         const auto rhFraction = hits_and_fractions[hitId].second;
0176         // if the fraction is zero or the hit does not belong to any SimTrackster,
0177         // set the SimTrackster Id for the hit to -1; this will
0178         // contribute to the number of noise hits
0179
0180         // MR Remove the case in which the fraction is 0, since this could be a
0181         // real hit that has been marked as halo.
0182         if (rhFraction == 0.) {
0183           hitsToSimTracksterId[hitId] = -2;
0184         }
0185         //Now check if there are SimTracksters linked to this rechit of the layercluster
0186         if (lc_find_in_ST == lcToSimTracksterId_Map.end()) {
0187           hitsToSimTracksterId[hitId] -= 1;
0188         } else {
0189           auto maxSTEnergyInLC = 0.f;
0190           auto maxSTId = -1;
0191           //Loop through all the linked SimTracksters
0192           for (const auto& st : lc_find_in_ST->second) {
0193             const auto stId = st.clusterId;
0194             STEnergyInLC[stId] += st.fraction * layerClusters[lcId].energy();
0195             // Keep track of which SimTrackster contributed the most, in terms
0196             // of energy, to this specific Layer Cluster.
0197             if (STEnergyInLC[stId] > maxSTEnergyInLC) {
0198               maxSTEnergyInLC = STEnergyInLC[stId];
0199               maxSTId = stId;
0200             }
0201           }
0202           hitsToSimTracksterId[hitId] = maxSTId;
0203         }
0204       }  // End loop over hits on a LayerCluster
0205
0206       for (const auto& c : hitsToSimTracksterId) {
0207         if (c < 0) {
0208           numberOfNoiseHitsInLC++;
0209         } else {
0210           occurrencesSTinLC[c]++;
0211         }
0212       }
0213
0214       for (const auto& c : occurrencesSTinLC) {
0215         if (c.second > maxSTNumberOfHitsInLC) {
0216           maxSTId_byNumberOfHits = c.first;
0217           maxSTNumberOfHitsInLC = c.second;
0218         }
0219       }
0220
0221       for (const auto& c : STEnergyInLC) {
0222         if (c.second > maxEnergySharedLCandST) {
0223           maxSTId_byEnergy = c.first;
0224           maxEnergySharedLCandST = c.second;
0225         }
0226       }
0227
0228       float totalSTEnergyOnLayer = 0.f;
0229       if (maxSTId_byEnergy >= 0) {
0230         totalSTEnergyOnLayer = tssInSimTrackster[maxSTId_byEnergy].energy;
0231         energyFractionOfSTinLC = maxEnergySharedLCandST / totalSTEnergyOnLayer;
0232         if (tracksters[tsId].raw_energy() > 0.f) {
0233           energyFractionOfLCinST = maxEnergySharedLCandST / tracksters[tsId].raw_energy();
0234         }
0235       }
0236
0237       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0238           << std::setw(12) << "TracksterID:\t" << std::setw(12) << "layerCluster\t" << std::setw(10) << "lc energy\t"
0239           << std::setw(5) << "nhits\t" << std::setw(12) << "noise hits\t" << std::setw(22) << "maxSTId_byNumberOfHits\t"
0240           << std::setw(8) << "nhitsST\t" << std::setw(13) << "maxSTId_byEnergy\t" << std::setw(20)
0241           << "maxEnergySharedLCandST\t" << std::setw(22) << "totalSTEnergyOnLayer\t" << std::setw(22)
0242           << "energyFractionOfLCinST\t" << std::setw(25) << "energyFractionOfSTinLC\t"
0243           << "\n";
0244       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0245           << std::setw(12) << tsId << "\t" << std::setw(12) << lcId << "\t" << std::setw(10)
0246           << tracksters[tsId].raw_energy() << "\t" << std::setw(5) << numberOfHitsInLC << "\t" << std::setw(12)
0247           << numberOfNoiseHitsInLC << "\t" << std::setw(22) << maxSTId_byNumberOfHits << "\t" << std::setw(8)
0248           << maxSTNumberOfHitsInLC << "\t" << std::setw(13) << maxSTId_byEnergy << "\t" << std::setw(20)
0249           << maxEnergySharedLCandST << "\t" << std::setw(22) << totalSTEnergyOnLayer << "\t" << std::setw(22)
0250           << energyFractionOfLCinST << "\t" << std::setw(25) << energyFractionOfSTinLC << "\n";
0251     }  // End of loop over LayerClusters in Trackster
0252   }  // End of loop over Tracksters
0253
0254   LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0255       << "Improved tssInSimTrackster INFO (Now containing the linked tracksters id and energy - score still empty)"
0256       << std::endl;
0257   for (size_t st = 0; st < tssInSimTrackster.size(); ++st) {
0258     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "For SimTrackster Idx: " << st << " we have: " << std::endl;
0259     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0260         << "    SimTracksterIdx: " << tssInSimTrackster[st].simTracksterId << std::endl;
0261     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "\tEnergy:\t" << tssInSimTrackster[st].energy << std::endl;
0262     double tot_energy = 0.;
0263     for (auto const& lcaf : tssInSimTrackster[st].lcs_and_fractions) {
0264       const auto lcEn = lcaf.second * layerClusters[lcaf.first].energy();
0265       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0266           << "\tLC/fraction/energy: " << (uint32_t)lcaf.first << "/" << lcaf.second << "/" << lcEn << std::endl;
0267       tot_energy += lcEn;
0268     }
0269     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "\tTot Sum lcaf: " << tot_energy << std::endl;
0270     for (auto const& ts : tssInSimTrackster[st].tracksterIdToEnergyAndScore) {
0271       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0272           << "\ttsIdx/energy/score: " << ts.first << "/" << ts.second.first << "/" << ts.second.second << std::endl;
0273     }
0274   }
0275
0276   LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "Improved lcToSimTracksterId_Map INFO" << std::endl;
0277   for (auto const& lc : lcToSimTracksterId_Map) {
0278     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0279         << "For lcId: " << (uint32_t)lc.first
0280         << " we have found the following connections with SimTracksters:" << std::endl;
0281     for (auto const& st : lc.second) {
0282       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0283           << "  SimTrackster Id: " << st.clusterId << " with fraction: " << st.fraction
0284           << " and energy: " << st.fraction * layerClusters[lc.first].energy() << std::endl;
0285     }
0286   }
0287 #endif
0288
0289   // Update stsInTrackster; compute the score Trackster-to-SimTrackster,
0290   // together with the returned AssociationMap
0291   for (unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
0292     // The SimTracksters contributing to the Trackster's LayerClusters should already be unique.
0293     // find the unique SimTracksters id contributing to the Trackster's LayerClusters
0294     std::sort(stsInTrackster[tsId].begin(), stsInTrackster[tsId].end());
0295     auto last = std::unique(stsInTrackster[tsId].begin(), stsInTrackster[tsId].end());
0296     stsInTrackster[tsId].erase(last, stsInTrackster[tsId].end());
0297
0298     // If a reconstructed Trackster has energy 0 but is linked to a
0299     // SimTrackster, assigned score 1
0300     if (tracksters[tsId].raw_energy() == 0. && !stsInTrackster[tsId].empty()) {
0301       for (auto& stPair : stsInTrackster[tsId]) {
0302         stPair.second.second = 1.;
0303         LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0304             << "TracksterId:\t " << tsId << "\tST id:\t" << stPair.first << "\tenergy" << stPair.second.first
0305             << "\tscore\t " << stPair.second.second << "\n";
0306       }
0307       continue;
0308     }
0309
0310     float invTracksterEnergyWeight = 0.f;
0311     for (unsigned int i = 0; i < tracksters[tsId].vertices().size(); ++i) {
0312       const auto lcId = tracksters[tsId].vertices(i);
0313       const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(i);
0314
0315       const auto& hits_and_fractions = layerClusters[lcId].hitsAndFractions();
0316       // Compute the correct normalization
0317       for (auto const& haf : hits_and_fractions) {
0318         const HGCRecHit* hit = hits_[hitMap_->at(haf.first)];
0319         invTracksterEnergyWeight += std::pow(lcFractionInTs * haf.second * hit->energy(), 2);
0320       }
0321     }
0322     invTracksterEnergyWeight = 1.f / invTracksterEnergyWeight;
0323
0324     for (unsigned int i = 0; i < tracksters[tsId].vertices().size(); ++i) {
0325       const auto lcId = tracksters[tsId].vertices(i);
0326       const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(i);
0327
0328       const bool lcWithST = (lcToSimTracksterId_Map.find(lcId) != lcToSimTracksterId_Map.end());
0329
0330       float lcEnergyWeight = pow(layerClusters[lcId].energy(), 2);
0331
0332       for (auto& stPair : stsInTrackster[tsId]) {
0333         float stFraction = 0.f;
0334         if (lcWithST) {
0335           const auto findLCIt = std::find(lcToSimTracksterId_Map[lcId].begin(),
0336                                           lcToSimTracksterId_Map[lcId].end(),
0337                                           ticl::lcInfoInTrackster{stPair.first, 0.f});
0338           if (findLCIt != lcToSimTracksterId_Map[lcId].end()) {
0339             stFraction = findLCIt->fraction;
0340           }
0341         }
0342         stPair.second.second +=
0343             (lcFractionInTs - stFraction) * (lcFractionInTs - stFraction) * lcEnergyWeight * invTracksterEnergyWeight;
0344 #ifdef EDM_ML_DEBUG
0345         LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0346             << "lcId:\t" << (uint32_t)lcId << "\ttracksterId:\t" << tsId << "\ttsFraction,stFraction:\t"
0347             << lcFractionInTs << ", " << stFraction << "\tlcEnergyWeight:\t" << lcEnergyWeight << "\tcurrent score:\t"
0348             << stPair.second.second << "\tinvTracksterEnergyWeight:\t" << invTracksterEnergyWeight << "\n";
0349 #endif
0350       }
0351     }  // End of loop over LayerClusters in Trackster
0352
0353 #ifdef EDM_ML_DEBUG
0354     if (stsInTrackster[tsId].empty())
0355       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "trackster Id:\t" << tsId << "\tST id:\t-1"
0356                                                        << "\tscore\t-1\n";
0357 #endif
0358   }  // End of loop over Tracksters
0359
0360   // Compute the SimTrackster-To-Trackster score
0361   for (const auto& stId : sTIndices) {
0362     float invSTEnergyWeight = 0.f;
0363
0364     const unsigned int STNumberOfLCs = tssInSimTrackster[stId].lcs_and_fractions.size();
0365     if (STNumberOfLCs == 0)
0366       continue;
0367
0368 #ifdef EDM_ML_DEBUG
0369     int tsWithMaxEnergyInST = -1;
0370     //energy of the most energetic TS from all that were linked to ST
0371     float maxEnergyTSinST = 0.f;
0372     float STenergy = tssInSimTrackster[stId].energy;
0373     //most energetic TS from all TSs linked to ST over ST energy.
0374     float STEnergyFractionInTS = 0.f;
0375     for (const auto& ts : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
0376       if (ts.second.first > maxEnergyTSinST) {
0377         maxEnergyTSinST = ts.second.first;
0378         tsWithMaxEnergyInST = ts.first;
0379       }
0380     }
0381     if (STenergy > 0.f)
0382       STEnergyFractionInTS = maxEnergyTSinST / STenergy;
0383
0384     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0385         << std::setw(12) << "simTrackster\t" << std::setw(15) << "st total energy\t" << std::setw(15)
0386         << "stEnergyOnLayer\t" << std::setw(14) << "STNhitsOnLayer\t" << std::setw(18) << "tsWithMaxEnergyInST\t"
0387         << std::setw(15) << "maxEnergyTSinST\t" << std::setw(20) << "STEnergyFractionInTS"
0388         << "\n";
0389     LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0390         << std::setw(12) << stId << "\t" << std::setw(15) << simTracksters[stId].raw_energy() << "\t" << std::setw(15)
0391         << STenergy << "\t" << std::setw(14) << STNumberOfLCs << "\t" << std::setw(18) << tsWithMaxEnergyInST << "\t"
0392         << std::setw(15) << maxEnergyTSinST << "\t" << std::setw(20) << STEnergyFractionInTS << "\n";
0393 #endif
0394
0395     // Compute the correct normalization
0396     for (auto const& lcaf : tssInSimTrackster[stId].lcs_and_fractions) {
0397       invSTEnergyWeight += std::pow(lcaf.second * layerClusters[lcaf.first].energy(), 2);
0398     }
0399     invSTEnergyWeight = 1.f / invSTEnergyWeight;
0400
0401     for (unsigned int i = 0; i < STNumberOfLCs; ++i) {
0402       auto& st_lcId = tssInSimTrackster[stId].lcs_and_fractions[i].first;
0403       auto& st_lcFraction = tssInSimTrackster[stId].lcs_and_fractions[i].second;
0404
0405       if (st_lcFraction == 0.f)
0406         continue;  // hopefully this should never happen
0407       float lcEnergyWeight = pow(layerClusters[st_lcId].energy(), 2);
0408       for (auto& tsPair : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
0409         unsigned int tsId = tsPair.first;
0410
0411         for (unsigned int i = 0; i < tracksters[tsId].vertices().size(); ++i) {
0412           const auto tsFraction = 1.f / tracksters[tsId].vertex_multiplicity(i);
0413
0414           tsPair.second.second +=
0415               (tsFraction - st_lcFraction) * (tsFraction - st_lcFraction) * lcEnergyWeight * invSTEnergyWeight;
0416 #ifdef EDM_ML_DEBUG
0417           LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0418               << "STLCId:\t" << (uint32_t)st_lcId << "\tTracksterId:\t" << tsId << "\t"
0419               << "tsFraction, stFraction:\t" << tsFraction << ", " << st_lcFraction << "\t"
0420               << "lcEnergyWeight:\t" << lcEnergyWeight << "\t"
0421               << "current score:\t" << tsPair.second.second << "\t"
0422               << "invSTEnergyWeight:\t" << invSTEnergyWeight << "\n";
0423 #endif
0424         }  // End of loop over Trackster's LayerClusters
0425       }  // End of loop over Tracksters linked to hits of this SimTrackster
0426     }  // End of loop over hits of SimTrackster on a Layer
0427 #ifdef EDM_ML_DEBUG
0428     if (tssInSimTrackster[stId].tracksterIdToEnergyAndScore.empty())
0429       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "ST Id:\t" << stId << "\tTS id:\t-1 "
0430                                                        << "\tscore\t-1\n";
0431
0432     for (const auto& tsPair : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
0433       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl")
0434           << "ST Id: \t" << stId << "\t TS id: \t" << tsPair.first << "\t score \t" << tsPair.second.second
0435           << "\t shared energy:\t" << tsPair.second.first << "\t shared energy fraction:\t"
0436           << (tsPair.second.first / STenergy) << "\n";
0437     }
0438 #endif
0439   }  // End loop over SimTrackster indices
0440   return {stsInTrackster, tssInSimTrackster};
0441 }
0442
0443 ticl::RecoToSimCollectionSimTracksters TSToSimTSAssociatorByEnergyScoreImpl::associateRecoToSim(
0444     const edm::Handle<ticl::TracksterCollection>& tCH,
0445     const edm::Handle<reco::CaloClusterCollection>& lCCH,
0446     const edm::Handle<ticl::TracksterCollection>& sTCH) const {
0447   ticl::RecoToSimCollectionSimTracksters returnValue(productGetter_);
0448   const auto& links = makeConnections(tCH, lCCH, sTCH);
0449
0450   const auto& stsInTrackster = std::get<0>(links);
0451   for (size_t tsId = 0; tsId < stsInTrackster.size(); ++tsId) {
0452     for (auto& stPair : stsInTrackster[tsId]) {
0453       LogDebug("TSToSimTSAssociatorByEnergyScoreImpl") << "Trackster Id:\t" << tsId << "\tSimTrackster id:\t"
0454                                                        << stPair.first << "\tscore:\t" << stPair.second.second << "\n";
0455       // Fill AssociationMap
0456       returnValue.insert(
0457           edm::Ref<ticl::TracksterCollection>(tCH, tsId),                            // Ref to TS
0458           std::make_pair(edm::Ref<ticl::TracksterCollection>(sTCH, stPair.first),    //Pair <Refo to TS>
0459                          std::make_pair(stPair.second.first, stPair.second.second))  // Pair <energy, score>
0460       );
0461     }
0462   }
0463   return returnValue;
0464 }
0465
0466 ticl::SimToRecoCollectionSimTracksters TSToSimTSAssociatorByEnergyScoreImpl::associateSimToReco(
0467     const edm::Handle<ticl::TracksterCollection>& tCH,
0468     const edm::Handle<reco::CaloClusterCollection>& lCCH,
0469     const edm::Handle<ticl::TracksterCollection>& sTCH) const {
0470   ticl::SimToRecoCollectionSimTracksters returnValue(productGetter_);
0471   const auto& links = makeConnections(tCH, lCCH, sTCH);
0472   const auto& tssInSimTrackster = std::get<1>(links);
0473   for (size_t stId = 0; stId < tssInSimTrackster.size(); ++stId) {
0474     for (auto& tsPair : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
0475       returnValue.insert(
0476           edm::Ref<ticl::TracksterCollection>(sTCH, stId),                           // Ref to ST
0477           std::make_pair(edm::Ref<ticl::TracksterCollection>(tCH, tsPair.first),     // Pair <Ref to TS,
0478                          std::make_pair(tsPair.second.first, tsPair.second.second))  // pair <energy, score> >
0479       );
0480     }
0481   }
0482   return returnValue;
0483 }