20 const auto& tracksters = *tCH.
product();
22 const auto& simClusters = *sCCH.
product();
23 auto nTracksters = tracksters.size();
27 auto nSimClusters = simClusters.size();
28 std::vector<size_t> sCIndices;
29 for (
unsigned int scId = 0; scId < nSimClusters; ++scId) {
31 simClusters[scId].g4Tracks()[0].eventId().bunchCrossing() != 0)) {
32 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
33 <<
"Excluding SimCluster from event: " << simClusters[scId].g4Tracks()[0].eventId().event()
34 <<
" with BX: " << simClusters[scId].g4Tracks()[0].eventId().bunchCrossing() << std::endl;
37 sCIndices.emplace_back(scId);
39 nSimClusters = sCIndices.size();
45 tssInSimCluster.resize(nSimClusters);
46 for (
unsigned int i = 0;
i < nSimClusters; ++
i) {
47 tssInSimCluster[
i].simClusterId =
i;
48 tssInSimCluster[
i].energy = 0.f;
49 tssInSimCluster[
i].hits_and_fractions.clear();
53 std::unordered_map<DetId, std::vector<hgcal::detIdInfoInCluster>> detIdToSimClusterId_Map;
54 for (
const auto& scId : sCIndices) {
55 const auto& hits_and_fractions = simClusters[scId].hits_and_fractions();
56 for (
const auto& it_haf : hits_and_fractions) {
57 const auto hitid = it_haf.first;
58 const auto itcheck =
hitMap_->find(hitid);
59 if (itcheck !=
hitMap_->end()) {
60 const auto hit_find_it = detIdToSimClusterId_Map.find(hitid);
61 if (hit_find_it == detIdToSimClusterId_Map.end()) {
62 detIdToSimClusterId_Map[hitid] = std::vector<hgcal::detIdInfoInCluster>();
64 detIdToSimClusterId_Map[hitid].emplace_back(scId, it_haf.second);
67 tssInSimCluster[scId].energy += it_haf.second *
hit->energy();
68 tssInSimCluster[scId].hits_and_fractions.emplace_back(hitid, it_haf.second);
74 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
75 <<
"tssInSimCluster INFO (Only SimCluster filled at the moment)" << std::endl;
76 for (
size_t sc = 0; sc < tssInSimCluster.size(); ++sc) {
77 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"For SimCluster Idx: " << sc <<
" we have: " << std::endl;
78 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
79 <<
"\tSimClusterIdx:\t" << tssInSimCluster[sc].simClusterId << std::endl;
80 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\tEnergy:\t" << tssInSimCluster[sc].energy << std::endl;
81 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\t# of clusters:\t" <<
layerClusters.size() << std::endl;
82 double tot_energy = 0.;
83 for (
auto const& haf : tssInSimCluster[sc].hits_and_fractions) {
84 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
85 <<
"\tHits/fraction/energy: " << (uint32_t)haf.first <<
"/" << haf.second <<
"/"
86 << haf.second *
hitMap_->at(haf.first)->energy() << std::endl;
87 tot_energy += haf.second *
hitMap_->at(haf.first)->energy();
89 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\tTot Sum haf: " << tot_energy << std::endl;
90 for (
auto const& ts : tssInSimCluster[sc].tracksterIdToEnergyAndScore) {
91 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
92 <<
"\ttsIdx/energy/score: " << ts.first <<
"/" << ts.second.first <<
"/" << ts.second.second << std::endl;
96 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"detIdToSimClusterId_Map INFO" << std::endl;
97 for (
auto const& detId : detIdToSimClusterId_Map) {
98 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
99 <<
"For detId: " << (uint32_t)detId.first
100 <<
" we have found the following connections with SimClusters:" << std::endl;
101 for (
auto const& sc : detId.second) {
102 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
103 <<
"\tSimCluster Id: " << sc.clusterId <<
" with fraction: " << sc.fraction
104 <<
" and energy: " << sc.fraction *
hitMap_->at(detId.first)->energy() << std::endl;
110 std::unordered_map<DetId, std::vector<hgcal::detIdInfoInCluster>> detIdToLayerClusterId_Map;
115 scsInTrackster.resize(nTracksters);
117 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
118 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
119 const auto lcId = tracksters[tsId].vertices(
i);
120 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
122 const std::vector<std::pair<DetId, float>>& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
123 unsigned int numberOfHitsInLC = hits_and_fractions.size();
125 for (
unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
126 const auto rh_detid = hits_and_fractions[hitId].first;
127 const auto rhFraction = hits_and_fractions[hitId].second;
129 const auto hit_find_in_LC = detIdToLayerClusterId_Map.find(rh_detid);
130 if (hit_find_in_LC == detIdToLayerClusterId_Map.end()) {
131 detIdToLayerClusterId_Map[rh_detid] = std::vector<hgcal::detIdInfoInCluster>();
133 detIdToLayerClusterId_Map[rh_detid].emplace_back(lcId, rhFraction);
135 const auto hit_find_in_SC = detIdToSimClusterId_Map.find(rh_detid);
137 if (hit_find_in_SC != detIdToSimClusterId_Map.end()) {
138 const auto itcheck =
hitMap_->find(rh_detid);
143 for (
const auto&
h : hit_find_in_SC->second) {
146 tssInSimCluster[
h.clusterId].tracksterIdToEnergyAndScore[tsId].first +=
147 lcFractionInTs *
h.fraction *
hit->energy();
149 scsInTrackster[tsId].emplace_back(
h.clusterId, 0.f);
157 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
158 for (
const auto& lcId : tracksters[tsId].
vertices()) {
159 const auto& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
160 unsigned int numberOfHitsInLC = hits_and_fractions.size();
170 std::vector<int> hitsToSimClusterId(numberOfHitsInLC);
172 int maxSCId_byNumberOfHits = -1;
174 unsigned int maxSCNumberOfHitsInLC = 0;
176 int maxSCId_byEnergy = -1;
178 float maxEnergySharedLCandSC = 0.f;
180 float energyFractionOfLCinSC = 0.f;
182 float energyFractionOfSCinLC = 0.f;
183 std::unordered_map<unsigned, unsigned> occurrencesSCinLC;
184 unsigned int numberOfNoiseHitsInLC = 0;
185 std::unordered_map<unsigned, float> SCEnergyInLC;
187 for (
unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
188 const auto rh_detid = hits_and_fractions[hitId].first;
189 const auto rhFraction = hits_and_fractions[hitId].second;
191 const auto hit_find_in_SC = detIdToSimClusterId_Map.find(rh_detid);
199 if (rhFraction == 0.) {
200 hitsToSimClusterId[hitId] = -2;
203 if (hit_find_in_SC == detIdToSimClusterId_Map.end()) {
204 hitsToSimClusterId[hitId] -= 1;
206 const auto itcheck =
hitMap_->find(rh_detid);
208 auto maxSCEnergyInLC = 0.f;
211 for (
const auto&
h : hit_find_in_SC->second) {
212 SCEnergyInLC[
h.clusterId] +=
h.fraction *
hit->energy();
215 if (SCEnergyInLC[
h.clusterId] > maxSCEnergyInLC) {
216 maxSCEnergyInLC = SCEnergyInLC[
h.clusterId];
217 maxSCId =
h.clusterId;
220 hitsToSimClusterId[hitId] = maxSCId;
224 for (
const auto&
c : hitsToSimClusterId) {
226 numberOfNoiseHitsInLC++;
228 occurrencesSCinLC[
c]++;
232 for (
const auto&
c : occurrencesSCinLC) {
233 if (
c.second > maxSCNumberOfHitsInLC) {
234 maxSCId_byNumberOfHits =
c.first;
235 maxSCNumberOfHitsInLC =
c.second;
239 for (
const auto&
c : SCEnergyInLC) {
240 if (
c.second > maxEnergySharedLCandSC) {
241 maxSCId_byEnergy =
c.first;
242 maxEnergySharedLCandSC =
c.second;
246 float totalSCEnergyOnLayer = 0.f;
247 if (maxSCId_byEnergy >= 0) {
248 totalSCEnergyOnLayer = tssInSimCluster[maxSCId_byEnergy].energy;
249 energyFractionOfSCinLC = maxEnergySharedLCandSC / totalSCEnergyOnLayer;
250 if (tracksters[tsId].raw_energy() > 0.
f) {
251 energyFractionOfLCinSC = maxEnergySharedLCandSC / tracksters[tsId].raw_energy();
255 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
256 << std::setw(12) <<
"TracksterID:\t" << std::setw(12) <<
"layerCluster\t" << std::setw(10) <<
"lc energy\t"
257 << std::setw(5) <<
"nhits\t" << std::setw(12) <<
"noise hits\t" << std::setw(22) <<
"maxSCId_byNumberOfHits\t"
258 << std::setw(8) <<
"nhitsSC\t" << std::setw(13) <<
"maxSCId_byEnergy\t" << std::setw(20)
259 <<
"maxEnergySharedLCandSC\t" << std::setw(22) <<
"totalSCEnergyOnLayer\t" << std::setw(22)
260 <<
"energyFractionOfLCinSC\t" << std::setw(25) <<
"energyFractionOfSCinLC\t"
262 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
263 << std::setw(12) << tsId <<
"\t" << std::setw(12) << lcId <<
"\t" << std::setw(10)
264 << tracksters[tsId].raw_energy() <<
"\t" << std::setw(5) << numberOfHitsInLC <<
"\t" << std::setw(12)
265 << numberOfNoiseHitsInLC <<
"\t" << std::setw(22) << maxSCId_byNumberOfHits <<
"\t" << std::setw(8)
266 << maxSCNumberOfHitsInLC <<
"\t" << std::setw(13) << maxSCId_byEnergy <<
"\t" << std::setw(20)
267 << maxEnergySharedLCandSC <<
"\t" << std::setw(22) << totalSCEnergyOnLayer <<
"\t" << std::setw(22)
268 << energyFractionOfLCinSC <<
"\t" << std::setw(25) << energyFractionOfSCinLC <<
"\n";
272 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
273 <<
"Improved tssInSimCluster INFO (Now containing the linked tracksters id and energy - score still empty)"
275 for (
size_t sc = 0; sc < tssInSimCluster.size(); ++sc) {
276 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"For SimCluster Idx: " << sc <<
" we have: " << std::endl;
277 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
278 <<
" SimClusterIdx: " << tssInSimCluster[sc].simClusterId << std::endl;
279 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\tEnergy:\t" << tssInSimCluster[sc].energy << std::endl;
280 double tot_energy = 0.;
281 for (
auto const& haf : tssInSimCluster[sc].hits_and_fractions) {
282 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
283 <<
"\tHits/fraction/energy: " << (uint32_t)haf.first <<
"/" << haf.second <<
"/"
284 << haf.second *
hitMap_->at(haf.first)->energy() << std::endl;
285 tot_energy += haf.second *
hitMap_->at(haf.first)->energy();
287 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\tTot Sum haf: " << tot_energy << std::endl;
288 for (
auto const& ts : tssInSimCluster[sc].tracksterIdToEnergyAndScore) {
289 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
290 <<
"\ttsIdx/energy/score: " << ts.first <<
"/" << ts.second.first <<
"/" << ts.second.second << std::endl;
294 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"Improved detIdToSimClusterId_Map INFO" << std::endl;
295 for (
auto const& sc : detIdToSimClusterId_Map) {
296 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
297 <<
"For detId: " << (uint32_t)sc.first
298 <<
" we have found the following connections with SimClusters:" << std::endl;
299 for (
auto const& sclu : sc.second) {
300 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
301 <<
" SimCluster Id: " << sclu.clusterId <<
" with fraction: " << sclu.fraction
302 <<
" and energy: " << sclu.fraction *
hitMap_->at(sc.first)->energy() << std::endl;
309 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
312 std::sort(scsInTrackster[tsId].begin(), scsInTrackster[tsId].
end());
314 scsInTrackster[tsId].erase(
last, scsInTrackster[tsId].
end());
318 if (tracksters[tsId].raw_energy() == 0. && !scsInTrackster[tsId].
empty()) {
319 for (
auto& scPair : scsInTrackster[tsId]) {
321 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
322 <<
"TracksterId:\t " << tsId <<
"\tSC id:\t" << scPair.first <<
"\tscore\t " << scPair.second <<
"\n";
327 float invTracksterEnergyWeight = 0.f;
328 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
329 const auto lcId = tracksters[tsId].vertices(
i);
330 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
332 const auto& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
334 for (
auto const& haf : hits_and_fractions) {
335 invTracksterEnergyWeight += (lcFractionInTs * haf.second *
hitMap_->at(haf.first)->energy()) *
336 (lcFractionInTs * haf.second *
hitMap_->at(haf.first)->energy());
339 invTracksterEnergyWeight = 1.f / invTracksterEnergyWeight;
341 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
342 const auto lcId = tracksters[tsId].vertices(
i);
343 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
345 const auto& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
346 unsigned int numberOfHitsInLC = hits_and_fractions.size();
347 for (
unsigned int i = 0;
i < numberOfHitsInLC; ++
i) {
348 DetId rh_detid = hits_and_fractions[
i].first;
349 float rhFraction = hits_and_fractions[
i].second * lcFractionInTs;
351 const bool hitWithSC = (detIdToSimClusterId_Map.find(rh_detid) != detIdToSimClusterId_Map.end());
353 const auto itcheck =
hitMap_->find(rh_detid);
355 float hitEnergyWeight =
hit->energy() *
hit->energy();
357 for (
auto& scPair : scsInTrackster[tsId]) {
358 float scFraction = 0.f;
360 const auto findHitIt =
std::find(detIdToSimClusterId_Map[rh_detid].begin(),
361 detIdToSimClusterId_Map[rh_detid].
end(),
363 if (findHitIt != detIdToSimClusterId_Map[rh_detid].
end())
364 scFraction = findHitIt->fraction;
367 (rhFraction - scFraction) * (rhFraction - scFraction) * hitEnergyWeight * invTracksterEnergyWeight;
369 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
370 <<
"rh_detid:\t" << (uint32_t)rh_detid <<
"\ttracksterId:\t" << tsId <<
"\t"
371 <<
"rhfraction,scFraction:\t" << rhFraction <<
", " << scFraction <<
"\t"
372 <<
"hitEnergyWeight:\t" << hitEnergyWeight <<
"\t"
373 <<
"current score:\t" << scPair.second <<
"\t"
374 <<
"invTracksterEnergyWeight:\t" << invTracksterEnergyWeight <<
"\n";
381 if (scsInTrackster[tsId].
empty())
382 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"trackster Id:\t" << tsId <<
"\tSC id:\t-1"
388 for (
const auto& scId : sCIndices) {
389 float invSCEnergyWeight = 0.f;
391 const unsigned int SCNumberOfHits = tssInSimCluster[scId].hits_and_fractions.size();
392 if (SCNumberOfHits == 0)
395 int tsWithMaxEnergyInSC = -1;
397 float maxEnergyTSinSC = 0.f;
398 float SCenergy = tssInSimCluster[scId].energy;
400 float SCEnergyFractionInTS = 0.f;
401 for (
const auto& ts : tssInSimCluster[scId].tracksterIdToEnergyAndScore) {
402 if (ts.second.first > maxEnergyTSinSC) {
403 maxEnergyTSinSC = ts.second.first;
404 tsWithMaxEnergyInSC = ts.first;
408 SCEnergyFractionInTS = maxEnergyTSinSC / SCenergy;
410 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
411 << std::setw(12) <<
"simcluster\t" << std::setw(15) <<
"sc total energy\t" << std::setw(15)
412 <<
"scEnergyOnLayer\t" << std::setw(14) <<
"SCNhitsOnLayer\t" << std::setw(18) <<
"tsWithMaxEnergyInSC\t"
413 << std::setw(15) <<
"maxEnergyTSinSC\t" << std::setw(20) <<
"SCEnergyFractionInTS"
415 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
416 << std::setw(12) << scId <<
"\t" << std::setw(15) << simClusters[scId].energy() <<
"\t" << std::setw(15)
417 << SCenergy <<
"\t" << std::setw(14) << SCNumberOfHits <<
"\t" << std::setw(18) << tsWithMaxEnergyInSC <<
"\t"
418 << std::setw(15) << maxEnergyTSinSC <<
"\t" << std::setw(20) << SCEnergyFractionInTS <<
"\n";
421 for (
auto const& haf : tssInSimCluster[scId].hits_and_fractions) {
422 invSCEnergyWeight +=
std::pow(haf.second *
hitMap_->at(haf.first)->energy(), 2);
424 invSCEnergyWeight = 1.f / invSCEnergyWeight;
426 for (
unsigned int i = 0;
i < SCNumberOfHits; ++
i) {
427 auto& sc_hitDetId = tssInSimCluster[scId].hits_and_fractions[
i].first;
428 auto& scFraction = tssInSimCluster[scId].hits_and_fractions[
i].second;
430 bool hitWithLC =
false;
431 if (scFraction == 0.
f)
433 const auto hit_find_in_LC = detIdToLayerClusterId_Map.find(sc_hitDetId);
434 if (hit_find_in_LC != detIdToLayerClusterId_Map.end())
436 const auto itcheck =
hitMap_->find(sc_hitDetId);
438 float hitEnergyWeight =
hit->energy() *
hit->energy();
439 for (
auto& tsPair : tssInSimCluster[scId].tracksterIdToEnergyAndScore) {
440 unsigned int tsId = tsPair.first;
441 float tsFraction = 0.f;
443 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
444 const auto lcId = tracksters[tsId].vertices(
i);
445 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
448 const auto findHitIt =
std::find(detIdToLayerClusterId_Map[sc_hitDetId].begin(),
449 detIdToLayerClusterId_Map[sc_hitDetId].
end(),
451 if (findHitIt != detIdToLayerClusterId_Map[sc_hitDetId].
end())
452 tsFraction = findHitIt->fraction * lcFractionInTs;
454 tsPair.second.second +=
455 (tsFraction - scFraction) * (tsFraction - scFraction) * hitEnergyWeight * invSCEnergyWeight;
457 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
458 <<
"SCDetId:\t" << (uint32_t)sc_hitDetId <<
"\tTracksterId:\t" << tsId <<
"\t"
459 <<
"tsFraction, scFraction:\t" << tsFraction <<
", " << scFraction <<
"\t"
460 <<
"hitEnergyWeight:\t" << hitEnergyWeight <<
"\t"
461 <<
"current score:\t" << tsPair.second.second <<
"\t"
462 <<
"invSCEnergyWeight:\t" << invSCEnergyWeight <<
"\n";
468 if (tssInSimCluster[scId].tracksterIdToEnergyAndScore.empty())
469 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"SC Id:\t" << scId <<
"\tTS id:\t-1 "
472 for (
const auto& tsPair : tssInSimCluster[scId].tracksterIdToEnergyAndScore) {
473 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
474 <<
"SC Id: \t" << scId <<
"\t TS id: \t" << tsPair.first <<
"\t score \t" << tsPair.second.second
475 <<
"\t shared energy:\t" << tsPair.second.first <<
"\t shared energy fraction:\t"
476 << (tsPair.second.first / SCenergy) <<
"\n";
480 return {scsInTrackster, tssInSimCluster};