19 const auto& tracksters = *tCH.
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);
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 +=
131 stsInTrackster[tsId].emplace_back(st.clusterId, std::make_pair(0.f, 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) {
291 stsInTrackster[tsId].erase(
last, stsInTrackster[tsId].
end());
295 if (tracksters[tsId].raw_energy() == 0. && !stsInTrackster[tsId].
empty()) {
296 for (
auto& stPair : stsInTrackster[tsId]) {
297 stPair.second.second = 1.;
298 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
299 <<
"TracksterId:\t " << tsId <<
"\tST id:\t" << stPair.first <<
"\tenergy" << stPair.second.first
300 <<
"\tscore\t " << stPair.second.second <<
"\n";
305 float invTracksterEnergyWeight = 0.f;
306 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
307 const auto lcId = tracksters[tsId].vertices(
i);
308 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
310 const auto& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
312 for (
auto const& haf : hits_and_fractions) {
313 invTracksterEnergyWeight +=
std::pow(lcFractionInTs * haf.second *
hitMap_->at(haf.first)->energy(), 2);
316 invTracksterEnergyWeight = 1.f / invTracksterEnergyWeight;
318 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
319 const auto lcId = tracksters[tsId].vertices(
i);
320 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
322 const bool lcWithST = (lcToSimTracksterId_Map.find(lcId) != lcToSimTracksterId_Map.end());
326 for (
auto& stPair : stsInTrackster[tsId]) {
327 float stFraction = 0.f;
329 const auto findLCIt =
std::find(lcToSimTracksterId_Map[lcId].
begin(),
330 lcToSimTracksterId_Map[lcId].
end(),
332 if (findLCIt != lcToSimTracksterId_Map[lcId].
end()) {
333 stFraction = findLCIt->fraction;
336 stPair.second.second +=
337 (lcFractionInTs - stFraction) * (lcFractionInTs - stFraction) * lcEnergyWeight * invTracksterEnergyWeight;
339 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
340 <<
"lcId:\t" << (uint32_t)lcId <<
"\ttracksterId:\t" << tsId <<
"\ttsFraction,stFraction:\t" 341 << lcFractionInTs <<
", " << stFraction <<
"\tlcEnergyWeight:\t" << lcEnergyWeight <<
"\tcurrent score:\t" 342 << stPair.second.second <<
"\tinvTracksterEnergyWeight:\t" << invTracksterEnergyWeight <<
"\n";
348 if (stsInTrackster[tsId].
empty())
349 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"trackster Id:\t" << tsId <<
"\tST id:\t-1" 355 for (
const auto& stId : sTIndices) {
356 float invSTEnergyWeight = 0.f;
358 const unsigned int STNumberOfLCs = tssInSimTrackster[stId].lcs_and_fractions.size();
359 if (STNumberOfLCs == 0)
363 int tsWithMaxEnergyInST = -1;
365 float maxEnergyTSinST = 0.f;
366 float STenergy = tssInSimTrackster[stId].energy;
368 float STEnergyFractionInTS = 0.f;
369 for (
const auto& ts : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
370 if (ts.second.first > maxEnergyTSinST) {
371 maxEnergyTSinST = ts.second.first;
372 tsWithMaxEnergyInST = ts.first;
376 STEnergyFractionInTS = maxEnergyTSinST / STenergy;
378 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
379 << std::setw(12) <<
"simTrackster\t" << std::setw(15) <<
"st total energy\t" << std::setw(15)
380 <<
"stEnergyOnLayer\t" << std::setw(14) <<
"STNhitsOnLayer\t" << std::setw(18) <<
"tsWithMaxEnergyInST\t" 381 << std::setw(15) <<
"maxEnergyTSinST\t" << std::setw(20) <<
"STEnergyFractionInTS" 383 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
384 << std::setw(12) << stId <<
"\t" << std::setw(15) << simTracksters[stId].raw_energy() <<
"\t" << std::setw(15)
385 << STenergy <<
"\t" << std::setw(14) << STNumberOfLCs <<
"\t" << std::setw(18) << tsWithMaxEnergyInST <<
"\t" 386 << std::setw(15) << maxEnergyTSinST <<
"\t" << std::setw(20) << STEnergyFractionInTS <<
"\n";
390 for (
auto const& lcaf : tssInSimTrackster[stId].lcs_and_fractions) {
393 invSTEnergyWeight = 1.f / invSTEnergyWeight;
395 for (
unsigned int i = 0;
i < STNumberOfLCs; ++
i) {
396 auto& st_lcId = tssInSimTrackster[stId].lcs_and_fractions[
i].first;
397 auto& st_lcFraction = tssInSimTrackster[stId].lcs_and_fractions[
i].second;
399 if (st_lcFraction == 0.
f)
402 for (
auto& tsPair : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
403 unsigned int tsId = tsPair.first;
405 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
406 const auto tsFraction = 1.f / tracksters[tsId].vertex_multiplicity(
i);
408 tsPair.second.second +=
409 (tsFraction - st_lcFraction) * (tsFraction - st_lcFraction) * lcEnergyWeight * invSTEnergyWeight;
411 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
412 <<
"STLCId:\t" << (uint32_t)st_lcId <<
"\tTracksterId:\t" << tsId <<
"\t" 413 <<
"tsFraction, stFraction:\t" << tsFraction <<
", " << st_lcFraction <<
"\t" 414 <<
"lcEnergyWeight:\t" << lcEnergyWeight <<
"\t" 415 <<
"current score:\t" << tsPair.second.second <<
"\t" 416 <<
"invSTEnergyWeight:\t" << invSTEnergyWeight <<
"\n";
422 if (tssInSimTrackster[stId].tracksterIdToEnergyAndScore.empty())
423 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"ST Id:\t" << stId <<
"\tTS id:\t-1 " 426 for (
const auto& tsPair : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
427 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
428 <<
"ST Id: \t" << stId <<
"\t TS id: \t" << tsPair.first <<
"\t score \t" << tsPair.second.second
429 <<
"\t shared energy:\t" << tsPair.second.first <<
"\t shared energy fraction:\t" 430 << (tsPair.second.first / STenergy) <<
"\n";
434 return {stsInTrackster, tssInSimTrackster};
for(int i=first, nt=offsets[nh];i< nt;i+=gridDim.x *blockDim.x)
T const * product() const
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
def unique(seq, keepstr=True)
std::vector< hgcal::simTracksterOnLayer > simTracksterToTrackster
std::vector< std::vector< std::pair< unsigned int, std::pair< float, float > > > > tracksterToSimTrackster
const std::unordered_map< DetId, const HGCRecHit * > * hitMap_
Power< A, B >::type pow(const A &a, const B &b)