19 const auto& tracksters = *tCH.
product();
20 const auto& layerClusters = *lCCH.
product();
21 const auto& simTracksters = *sTCH.
product();
22 auto nTracksters = tracksters.size();
25 auto nSimTracksters = simTracksters.size();
26 std::vector<size_t> sTIndices;
27 for (
unsigned int stId = 0; stId < nSimTracksters; ++stId) {
29 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
30 <<
"Excluding SimTrackster " << stId <<
" witH no vertices!" << std::endl;
33 sTIndices.emplace_back(stId);
35 nSimTracksters = sTIndices.size();
42 tssInSimTrackster.resize(nSimTracksters);
43 for (
unsigned int i = 0;
i < nSimTracksters; ++
i) {
44 tssInSimTrackster[
i].simTracksterId =
i;
45 tssInSimTrackster[
i].energy = 0.f;
46 tssInSimTrackster[
i].lcs_and_fractions.clear();
53 std::unordered_map<int, std::vector<hgcal::lcInfoInTrackster>> lcToSimTracksterId_Map;
54 for (
const auto& stId : sTIndices) {
55 const auto& lcs = simTracksters[stId].vertices();
57 for (
const auto& lcId : lcs) {
58 const auto fraction = 1.f / simTracksters[stId].vertex_multiplicity(lcInSimTst++);
60 const auto lc_find_it = lcToSimTracksterId_Map.find(lcId);
61 if (lc_find_it == lcToSimTracksterId_Map.end()) {
62 lcToSimTracksterId_Map[lcId] = std::vector<hgcal::lcInfoInTrackster>();
64 lcToSimTracksterId_Map[lcId].emplace_back(stId,
fraction);
66 tssInSimTrackster[stId].energy +=
fraction * layerClusters[lcId].energy();
67 tssInSimTrackster[stId].lcs_and_fractions.emplace_back(lcId,
fraction);
72 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
73 <<
"tssInSimTrackster INFO (Only SimTrackster filled at the moment)" << std::endl;
74 for (
size_t st = 0; st < tssInSimTrackster.size(); ++st) {
75 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"For SimTrackster Idx: " << st <<
" we have: " << std::endl;
76 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
77 <<
"\tSimTracksterIdx:\t" << tssInSimTrackster[st].simTracksterId << std::endl;
78 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"\tEnergy:\t" << tssInSimTrackster[st].energy << std::endl;
79 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"\t# of clusters:\t" << layerClusters.size() << std::endl;
80 double tot_energy = 0.;
81 for (
auto const& lcaf : tssInSimTrackster[st].lcs_and_fractions) {
82 const auto lcEn = lcaf.second * layerClusters[lcaf.first].energy();
83 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
84 <<
"\tLC/fraction/energy: " << (uint32_t)lcaf.first <<
"/" << lcaf.second <<
"/" << lcEn << std::endl;
87 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"\tTot Sum lcaf: " << tot_energy << std::endl;
88 for (
auto const& ts : tssInSimTrackster[st].tracksterIdToEnergyAndScore) {
89 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
90 <<
"\ttsIdx/energy/score: " << ts.first <<
"/" << ts.second.first <<
"/" << ts.second.second << std::endl;
94 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"lcToSimTracksterId_Map INFO" << std::endl;
95 for (
auto const& lc : lcToSimTracksterId_Map) {
96 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
97 <<
"For lcId: " << (uint32_t)lc.first
98 <<
" we have found the following connections with SimTracksters:" << std::endl;
99 for (
auto const& st : lc.second) {
100 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
101 <<
"\tSimTrackster Id: " << st.clusterId <<
" with fraction: " << st.fraction
102 <<
" and energy: " << st.fraction * layerClusters[lc.first].energy() << std::endl;
112 stsInTrackster.resize(nTracksters);
114 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
115 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
116 const auto lcId = tracksters[tsId].vertices(
i);
117 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
119 const auto lc_find_in_ST = lcToSimTracksterId_Map.find(lcId);
121 if (lc_find_in_ST != lcToSimTracksterId_Map.end()) {
125 for (
const auto& st : lc_find_in_ST->second) {
128 tssInSimTrackster[st.clusterId].tracksterIdToEnergyAndScore[tsId].first +=
129 lcFractionInTs * st.fraction * layerClusters[lcId].energy();
131 stsInTrackster[tsId].emplace_back(st.clusterId, 0.f);
138 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
139 for (
const auto& lcId : tracksters[tsId].
vertices()) {
140 const auto& hits_and_fractions = layerClusters[lcId].hitsAndFractions();
141 unsigned int numberOfHitsInLC = hits_and_fractions.size();
151 std::vector<int> hitsToSimTracksterId(numberOfHitsInLC);
153 int maxSTId_byNumberOfHits = -1;
155 unsigned int maxSTNumberOfHitsInLC = 0;
157 int maxSTId_byEnergy = -1;
159 float maxEnergySharedLCandST = 0.f;
161 float energyFractionOfLCinST = 0.f;
163 float energyFractionOfSTinLC = 0.f;
164 std::unordered_map<unsigned, unsigned> occurrencesSTinLC;
165 unsigned int numberOfNoiseHitsInLC = 0;
166 std::unordered_map<unsigned, float> STEnergyInLC;
168 const auto lc_find_in_ST = lcToSimTracksterId_Map.find(lcId);
169 for (
unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
170 const auto rhFraction = hits_and_fractions[hitId].second;
177 if (rhFraction == 0.) {
178 hitsToSimTracksterId[hitId] = -2;
181 if (lc_find_in_ST == lcToSimTracksterId_Map.end()) {
182 hitsToSimTracksterId[hitId] -= 1;
184 auto maxSTEnergyInLC = 0.f;
187 for (
const auto& st : lc_find_in_ST->second) {
188 const auto stId = st.clusterId;
189 STEnergyInLC[stId] += st.fraction * layerClusters[lcId].energy();
192 if (STEnergyInLC[stId] > maxSTEnergyInLC) {
193 maxSTEnergyInLC = STEnergyInLC[stId];
197 hitsToSimTracksterId[hitId] = maxSTId;
201 for (
const auto&
c : hitsToSimTracksterId) {
203 numberOfNoiseHitsInLC++;
205 occurrencesSTinLC[
c]++;
209 for (
const auto&
c : occurrencesSTinLC) {
210 if (
c.second > maxSTNumberOfHitsInLC) {
211 maxSTId_byNumberOfHits =
c.first;
212 maxSTNumberOfHitsInLC =
c.second;
216 for (
const auto&
c : STEnergyInLC) {
217 if (
c.second > maxEnergySharedLCandST) {
218 maxSTId_byEnergy =
c.first;
219 maxEnergySharedLCandST =
c.second;
223 float totalSTEnergyOnLayer = 0.f;
224 if (maxSTId_byEnergy >= 0) {
225 totalSTEnergyOnLayer = tssInSimTrackster[maxSTId_byEnergy].energy;
226 energyFractionOfSTinLC = maxEnergySharedLCandST / totalSTEnergyOnLayer;
227 if (tracksters[tsId].raw_energy() > 0.
f) {
228 energyFractionOfLCinST = maxEnergySharedLCandST / tracksters[tsId].raw_energy();
232 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
233 << std::setw(12) <<
"TracksterID:\t" << std::setw(12) <<
"layerCluster\t" << std::setw(10) <<
"lc energy\t"
234 << std::setw(5) <<
"nhits\t" << std::setw(12) <<
"noise hits\t" << std::setw(22) <<
"maxSTId_byNumberOfHits\t"
235 << std::setw(8) <<
"nhitsST\t" << std::setw(13) <<
"maxSTId_byEnergy\t" << std::setw(20)
236 <<
"maxEnergySharedLCandST\t" << std::setw(22) <<
"totalSTEnergyOnLayer\t" << std::setw(22)
237 <<
"energyFractionOfLCinST\t" << std::setw(25) <<
"energyFractionOfSTinLC\t"
239 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
240 << std::setw(12) << tsId <<
"\t" << std::setw(12) << lcId <<
"\t" << std::setw(10)
241 << tracksters[tsId].raw_energy() <<
"\t" << std::setw(5) << numberOfHitsInLC <<
"\t" << std::setw(12)
242 << numberOfNoiseHitsInLC <<
"\t" << std::setw(22) << maxSTId_byNumberOfHits <<
"\t" << std::setw(8)
243 << maxSTNumberOfHitsInLC <<
"\t" << std::setw(13) << maxSTId_byEnergy <<
"\t" << std::setw(20)
244 << maxEnergySharedLCandST <<
"\t" << std::setw(22) << totalSTEnergyOnLayer <<
"\t" << std::setw(22)
245 << energyFractionOfLCinST <<
"\t" << std::setw(25) << energyFractionOfSTinLC <<
"\n";
249 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
250 <<
"Improved tssInSimTrackster INFO (Now containing the linked tracksters id and energy - score still empty)"
252 for (
size_t st = 0; st < tssInSimTrackster.size(); ++st) {
253 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"For SimTrackster Idx: " << st <<
" we have: " << std::endl;
254 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
255 <<
" SimTracksterIdx: " << tssInSimTrackster[st].simTracksterId << std::endl;
256 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"\tEnergy:\t" << tssInSimTrackster[st].energy << std::endl;
257 double tot_energy = 0.;
258 for (
auto const& lcaf : tssInSimTrackster[st].lcs_and_fractions) {
259 const auto lcEn = lcaf.second * layerClusters[lcaf.first].energy();
260 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
261 <<
"\tLC/fraction/energy: " << (uint32_t)lcaf.first <<
"/" << lcaf.second <<
"/" << lcEn << std::endl;
264 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"\tTot Sum lcaf: " << tot_energy << std::endl;
265 for (
auto const& ts : tssInSimTrackster[st].tracksterIdToEnergyAndScore) {
266 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
267 <<
"\ttsIdx/energy/score: " << ts.first <<
"/" << ts.second.first <<
"/" << ts.second.second << std::endl;
271 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"Improved lcToSimTracksterId_Map INFO" << std::endl;
272 for (
auto const& lc : lcToSimTracksterId_Map) {
273 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
274 <<
"For lcId: " << (uint32_t)lc.first
275 <<
" we have found the following connections with SimTracksters:" << std::endl;
276 for (
auto const& st : lc.second) {
277 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
278 <<
" SimTrackster Id: " << st.clusterId <<
" with fraction: " << st.fraction
279 <<
" and energy: " << st.fraction * layerClusters[lc.first].energy() << std::endl;
286 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
289 std::sort(stsInTrackster[tsId].
begin(), stsInTrackster[tsId].
end());
291 stsInTrackster[tsId].erase(
last, stsInTrackster[tsId].
end());
295 if (tracksters[tsId].raw_energy() == 0. && !stsInTrackster[tsId].
empty()) {
296 for (
auto& stPair : stsInTrackster[tsId]) {
298 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
299 <<
"TracksterId:\t " << tsId <<
"\tST id:\t" << stPair.first <<
"\tscore\t " << stPair.second <<
"\n";
304 float invTracksterEnergyWeight = 0.f;
305 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
306 const auto lcId = tracksters[tsId].vertices(
i);
307 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
309 const auto& hits_and_fractions = layerClusters[lcId].hitsAndFractions();
311 for (
auto const& haf : hits_and_fractions) {
312 invTracksterEnergyWeight +=
std::pow(lcFractionInTs * haf.second *
hitMap_->at(haf.first)->energy(), 2);
315 invTracksterEnergyWeight = 1.f / invTracksterEnergyWeight;
317 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
318 const auto lcId = tracksters[tsId].vertices(
i);
319 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
321 const bool lcWithST = (lcToSimTracksterId_Map.find(lcId) != lcToSimTracksterId_Map.end());
323 float lcEnergyWeight =
pow(layerClusters[lcId].
energy(), 2);
325 for (
auto& stPair : stsInTrackster[tsId]) {
326 float stFraction = 0.f;
328 const auto findLCIt =
std::find(lcToSimTracksterId_Map[lcId].
begin(),
329 lcToSimTracksterId_Map[lcId].
end(),
331 if (findLCIt != lcToSimTracksterId_Map[lcId].
end()) {
332 stFraction = findLCIt->fraction;
336 (lcFractionInTs - stFraction) * (lcFractionInTs - stFraction) * lcEnergyWeight * invTracksterEnergyWeight;
338 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
339 <<
"lcId:\t" << (uint32_t)lcId <<
"\ttracksterId:\t" << tsId <<
"\ttsFraction,stFraction:\t"
340 << lcFractionInTs <<
", " << stFraction <<
"\tlcEnergyWeight:\t" << lcEnergyWeight <<
"\tcurrent score:\t"
341 << stPair.second <<
"\tinvTracksterEnergyWeight:\t" << invTracksterEnergyWeight <<
"\n";
347 if (stsInTrackster[tsId].
empty())
348 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"trackster Id:\t" << tsId <<
"\tST id:\t-1"
354 for (
const auto& stId : sTIndices) {
355 float invSTEnergyWeight = 0.f;
357 const unsigned int STNumberOfLCs = tssInSimTrackster[stId].lcs_and_fractions.size();
358 if (STNumberOfLCs == 0)
362 int tsWithMaxEnergyInST = -1;
364 float maxEnergyTSinST = 0.f;
365 float STenergy = tssInSimTrackster[stId].energy;
367 float STEnergyFractionInTS = 0.f;
368 for (
const auto& ts : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
369 if (ts.second.first > maxEnergyTSinST) {
370 maxEnergyTSinST = ts.second.first;
371 tsWithMaxEnergyInST = ts.first;
375 STEnergyFractionInTS = maxEnergyTSinST / STenergy;
377 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
378 << std::setw(12) <<
"simTrackster\t" << std::setw(15) <<
"st total energy\t" << std::setw(15)
379 <<
"stEnergyOnLayer\t" << std::setw(14) <<
"STNhitsOnLayer\t" << std::setw(18) <<
"tsWithMaxEnergyInST\t"
380 << std::setw(15) <<
"maxEnergyTSinST\t" << std::setw(20) <<
"STEnergyFractionInTS"
382 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
383 << std::setw(12) << stId <<
"\t" << std::setw(15) << simTracksters[stId].raw_energy() <<
"\t" << std::setw(15)
384 << STenergy <<
"\t" << std::setw(14) << STNumberOfLCs <<
"\t" << std::setw(18) << tsWithMaxEnergyInST <<
"\t"
385 << std::setw(15) << maxEnergyTSinST <<
"\t" << std::setw(20) << STEnergyFractionInTS <<
"\n";
389 for (
auto const& lcaf : tssInSimTrackster[stId].lcs_and_fractions) {
390 invSTEnergyWeight +=
std::pow(lcaf.second * layerClusters[lcaf.first].energy(), 2);
392 invSTEnergyWeight = 1.f / invSTEnergyWeight;
394 for (
unsigned int i = 0;
i < STNumberOfLCs; ++
i) {
395 auto& st_lcId = tssInSimTrackster[stId].lcs_and_fractions[
i].first;
396 auto& st_lcFraction = tssInSimTrackster[stId].lcs_and_fractions[
i].second;
398 if (st_lcFraction == 0.
f)
400 float lcEnergyWeight =
pow(layerClusters[st_lcId].
energy(), 2);
401 for (
auto& tsPair : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
402 unsigned int tsId = tsPair.first;
404 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
405 const auto tsFraction = 1.f / tracksters[tsId].vertex_multiplicity(
i);
407 tsPair.second.second +=
408 (tsFraction - st_lcFraction) * (tsFraction - st_lcFraction) * lcEnergyWeight * invSTEnergyWeight;
410 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
411 <<
"STLCId:\t" << (uint32_t)st_lcId <<
"\tTracksterId:\t" << tsId <<
"\t"
412 <<
"tsFraction, stFraction:\t" << tsFraction <<
", " << st_lcFraction <<
"\t"
413 <<
"lcEnergyWeight:\t" << lcEnergyWeight <<
"\t"
414 <<
"current score:\t" << tsPair.second.second <<
"\t"
415 <<
"invSTEnergyWeight:\t" << invSTEnergyWeight <<
"\n";
421 if (tssInSimTrackster[stId].tracksterIdToEnergyAndScore.empty())
422 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"ST Id:\t" << stId <<
"\tTS id:\t-1 "
425 for (
const auto& tsPair : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
426 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
427 <<
"ST Id: \t" << stId <<
"\t TS id: \t" << tsPair.first <<
"\t score \t" << tsPair.second.second
428 <<
"\t shared energy:\t" << tsPair.second.first <<
"\t shared energy fraction:\t"
429 << (tsPair.second.first / STenergy) <<
"\n";
433 return {stsInTrackster, tssInSimTrackster};
const edm::EventSetup & c
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
std::vector< std::vector< std::pair< unsigned int, float > > > tracksterToSimTrackster
std::vector< hgcal::simTracksterOnLayer > simTracksterToTrackster
for(Iditer=Id.begin();Iditer!=Id.end();Iditer++)
T const * product() const
const std::unordered_map< DetId, const HGCRecHit * > * hitMap_
Power< A, B >::type pow(const A &a, const B &b)