24 const auto& tracksters = *tCH.
product();
26 const auto& simTracksters = *sTCH.
product();
27 auto nTracksters = tracksters.size();
30 auto nSimTracksters = simTracksters.size();
31 std::vector<size_t> sTIndices;
32 for (
unsigned int stId = 0; stId < nSimTracksters; ++stId) {
34 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
35 <<
"Excluding SimTrackster " << stId <<
" witH no vertices!" << std::endl;
38 sTIndices.emplace_back(stId);
40 nSimTracksters = sTIndices.size();
47 tssInSimTrackster.resize(nSimTracksters);
48 for (
unsigned int i = 0;
i < nSimTracksters; ++
i) {
49 tssInSimTrackster[
i].simTracksterId =
i;
50 tssInSimTrackster[
i].energy = 0.f;
51 tssInSimTrackster[
i].lcs_and_fractions.clear();
58 std::unordered_map<int, std::vector<ticl::lcInfoInTrackster>> lcToSimTracksterId_Map;
59 for (
const auto& stId : sTIndices) {
60 const auto& lcs = simTracksters[stId].vertices();
62 for (
const auto& lcId : lcs) {
63 const auto fraction = 1.f / simTracksters[stId].vertex_multiplicity(lcInSimTst++);
65 const auto lc_find_it = lcToSimTracksterId_Map.find(lcId);
66 if (lc_find_it == lcToSimTracksterId_Map.end()) {
67 lcToSimTracksterId_Map[lcId] = std::vector<ticl::lcInfoInTrackster>();
69 lcToSimTracksterId_Map[lcId].emplace_back(stId,
fraction);
72 tssInSimTrackster[stId].lcs_and_fractions.emplace_back(lcId,
fraction);
77 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
78 <<
"tssInSimTrackster INFO (Only SimTrackster filled at the moment)" << std::endl;
79 for (
size_t st = 0; st < tssInSimTrackster.size(); ++st) {
80 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"For SimTrackster Idx: " << st <<
" we have: " << std::endl;
81 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
82 <<
"\tSimTracksterIdx:\t" << tssInSimTrackster[st].simTracksterId << std::endl;
83 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"\tEnergy:\t" << tssInSimTrackster[st].energy << std::endl;
84 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"\t# of clusters:\t" <<
layerClusters.size() << std::endl;
85 double tot_energy = 0.;
86 for (
auto const& lcaf : tssInSimTrackster[st].lcs_and_fractions) {
87 const auto lcEn = lcaf.second *
layerClusters[lcaf.first].energy();
88 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
89 <<
"\tLC/fraction/energy: " << (uint32_t)lcaf.first <<
"/" << lcaf.second <<
"/" << lcEn << std::endl;
92 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"\tTot Sum lcaf: " << tot_energy << std::endl;
93 for (
auto const& ts : tssInSimTrackster[st].tracksterIdToEnergyAndScore) {
94 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
95 <<
"\ttsIdx/energy/score: " << ts.first <<
"/" << ts.second.first <<
"/" << ts.second.second << std::endl;
99 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"lcToSimTracksterId_Map INFO" << std::endl;
100 for (
auto const& lc : lcToSimTracksterId_Map) {
101 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
102 <<
"For lcId: " << (uint32_t)lc.first
103 <<
" we have found the following connections with SimTracksters:" << std::endl;
104 for (
auto const& st : lc.second) {
105 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
106 <<
"\tSimTrackster Id: " << st.clusterId <<
" with fraction: " << st.fraction
107 <<
" and energy: " << st.fraction *
layerClusters[lc.first].energy() << std::endl;
117 stsInTrackster.resize(nTracksters);
119 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
120 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
121 const auto lcId = tracksters[tsId].vertices(
i);
122 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
124 const auto lc_find_in_ST = lcToSimTracksterId_Map.find(lcId);
126 if (lc_find_in_ST != lcToSimTracksterId_Map.end()) {
130 for (
const auto& st : lc_find_in_ST->second) {
133 tssInSimTrackster[st.clusterId].tracksterIdToEnergyAndScore[tsId].first +=
136 stsInTrackster[tsId].emplace_back(st.clusterId, std::make_pair(0.f, 0.f));
143 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
144 for (
const auto& lcId : tracksters[tsId].
vertices()) {
145 const auto& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
146 unsigned int numberOfHitsInLC = hits_and_fractions.size();
156 std::vector<int> hitsToSimTracksterId(numberOfHitsInLC);
158 int maxSTId_byNumberOfHits = -1;
160 unsigned int maxSTNumberOfHitsInLC = 0;
162 int maxSTId_byEnergy = -1;
164 float maxEnergySharedLCandST = 0.f;
166 float energyFractionOfLCinST = 0.f;
168 float energyFractionOfSTinLC = 0.f;
169 std::unordered_map<unsigned, unsigned> occurrencesSTinLC;
170 unsigned int numberOfNoiseHitsInLC = 0;
171 std::unordered_map<unsigned, float> STEnergyInLC;
173 const auto lc_find_in_ST = lcToSimTracksterId_Map.find(lcId);
174 for (
unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
175 const auto rhFraction = hits_and_fractions[hitId].second;
182 if (rhFraction == 0.) {
183 hitsToSimTracksterId[hitId] = -2;
186 if (lc_find_in_ST == lcToSimTracksterId_Map.end()) {
187 hitsToSimTracksterId[hitId] -= 1;
189 auto maxSTEnergyInLC = 0.f;
192 for (
const auto& st : lc_find_in_ST->second) {
193 const auto stId = st.clusterId;
194 STEnergyInLC[stId] += st.fraction *
layerClusters[lcId].energy();
197 if (STEnergyInLC[stId] > maxSTEnergyInLC) {
198 maxSTEnergyInLC = STEnergyInLC[stId];
202 hitsToSimTracksterId[hitId] = maxSTId;
206 for (
const auto&
c : hitsToSimTracksterId) {
208 numberOfNoiseHitsInLC++;
210 occurrencesSTinLC[
c]++;
214 for (
const auto&
c : occurrencesSTinLC) {
215 if (
c.second > maxSTNumberOfHitsInLC) {
216 maxSTId_byNumberOfHits =
c.first;
217 maxSTNumberOfHitsInLC =
c.second;
221 for (
const auto&
c : STEnergyInLC) {
222 if (
c.second > maxEnergySharedLCandST) {
223 maxSTId_byEnergy =
c.first;
224 maxEnergySharedLCandST =
c.second;
228 float totalSTEnergyOnLayer = 0.f;
229 if (maxSTId_byEnergy >= 0) {
230 totalSTEnergyOnLayer = tssInSimTrackster[maxSTId_byEnergy].energy;
231 energyFractionOfSTinLC = maxEnergySharedLCandST / totalSTEnergyOnLayer;
232 if (tracksters[tsId].raw_energy() > 0.
f) {
233 energyFractionOfLCinST = maxEnergySharedLCandST / tracksters[tsId].raw_energy();
237 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
238 << std::setw(12) <<
"TracksterID:\t" << std::setw(12) <<
"layerCluster\t" << std::setw(10) <<
"lc energy\t" 239 << std::setw(5) <<
"nhits\t" << std::setw(12) <<
"noise hits\t" << std::setw(22) <<
"maxSTId_byNumberOfHits\t" 240 << std::setw(8) <<
"nhitsST\t" << std::setw(13) <<
"maxSTId_byEnergy\t" << std::setw(20)
241 <<
"maxEnergySharedLCandST\t" << std::setw(22) <<
"totalSTEnergyOnLayer\t" << std::setw(22)
242 <<
"energyFractionOfLCinST\t" << std::setw(25) <<
"energyFractionOfSTinLC\t" 244 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
245 << std::setw(12) << tsId <<
"\t" << std::setw(12) << lcId <<
"\t" << std::setw(10)
246 << tracksters[tsId].raw_energy() <<
"\t" << std::setw(5) << numberOfHitsInLC <<
"\t" << std::setw(12)
247 << numberOfNoiseHitsInLC <<
"\t" << std::setw(22) << maxSTId_byNumberOfHits <<
"\t" << std::setw(8)
248 << maxSTNumberOfHitsInLC <<
"\t" << std::setw(13) << maxSTId_byEnergy <<
"\t" << std::setw(20)
249 << maxEnergySharedLCandST <<
"\t" << std::setw(22) << totalSTEnergyOnLayer <<
"\t" << std::setw(22)
250 << energyFractionOfLCinST <<
"\t" << std::setw(25) << energyFractionOfSTinLC <<
"\n";
254 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
255 <<
"Improved tssInSimTrackster INFO (Now containing the linked tracksters id and energy - score still empty)" 257 for (
size_t st = 0; st < tssInSimTrackster.size(); ++st) {
258 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"For SimTrackster Idx: " << st <<
" we have: " << std::endl;
259 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
260 <<
" SimTracksterIdx: " << tssInSimTrackster[st].simTracksterId << std::endl;
261 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"\tEnergy:\t" << tssInSimTrackster[st].energy << std::endl;
262 double tot_energy = 0.;
263 for (
auto const& lcaf : tssInSimTrackster[st].lcs_and_fractions) {
264 const auto lcEn = lcaf.second *
layerClusters[lcaf.first].energy();
265 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
266 <<
"\tLC/fraction/energy: " << (uint32_t)lcaf.first <<
"/" << lcaf.second <<
"/" << lcEn << std::endl;
269 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"\tTot Sum lcaf: " << tot_energy << std::endl;
270 for (
auto const& ts : tssInSimTrackster[st].tracksterIdToEnergyAndScore) {
271 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
272 <<
"\ttsIdx/energy/score: " << ts.first <<
"/" << ts.second.first <<
"/" << ts.second.second << std::endl;
276 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"Improved lcToSimTracksterId_Map INFO" << std::endl;
277 for (
auto const& lc : lcToSimTracksterId_Map) {
278 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
279 <<
"For lcId: " << (uint32_t)lc.first
280 <<
" we have found the following connections with SimTracksters:" << std::endl;
281 for (
auto const& st : lc.second) {
282 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
283 <<
" SimTrackster Id: " << st.clusterId <<
" with fraction: " << st.fraction
284 <<
" and energy: " << st.fraction *
layerClusters[lc.first].energy() << std::endl;
291 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
296 stsInTrackster[tsId].erase(
last, stsInTrackster[tsId].
end());
300 if (tracksters[tsId].raw_energy() == 0. && !stsInTrackster[tsId].
empty()) {
301 for (
auto& stPair : stsInTrackster[tsId]) {
302 stPair.second.second = 1.;
303 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
304 <<
"TracksterId:\t " << tsId <<
"\tST id:\t" << stPair.first <<
"\tenergy" << stPair.second.first
305 <<
"\tscore\t " << stPair.second.second <<
"\n";
310 float invTracksterEnergyWeight = 0.f;
311 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
312 const auto lcId = tracksters[tsId].vertices(
i);
313 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
315 const auto& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
317 for (
auto const& haf : hits_and_fractions) {
319 invTracksterEnergyWeight +=
std::pow(lcFractionInTs * haf.second *
hit->energy(), 2);
322 invTracksterEnergyWeight = 1.f / invTracksterEnergyWeight;
324 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
325 const auto lcId = tracksters[tsId].vertices(
i);
326 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
328 const bool lcWithST = (lcToSimTracksterId_Map.find(lcId) != lcToSimTracksterId_Map.end());
332 for (
auto& stPair : stsInTrackster[tsId]) {
333 float stFraction = 0.f;
335 const auto findLCIt =
std::find(lcToSimTracksterId_Map[lcId].
begin(),
336 lcToSimTracksterId_Map[lcId].
end(),
338 if (findLCIt != lcToSimTracksterId_Map[lcId].
end()) {
339 stFraction = findLCIt->fraction;
342 stPair.second.second +=
343 (lcFractionInTs - stFraction) * (lcFractionInTs - stFraction) * lcEnergyWeight * invTracksterEnergyWeight;
345 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
346 <<
"lcId:\t" << (uint32_t)lcId <<
"\ttracksterId:\t" << tsId <<
"\ttsFraction,stFraction:\t" 347 << lcFractionInTs <<
", " << stFraction <<
"\tlcEnergyWeight:\t" << lcEnergyWeight <<
"\tcurrent score:\t" 348 << stPair.second.second <<
"\tinvTracksterEnergyWeight:\t" << invTracksterEnergyWeight <<
"\n";
354 if (stsInTrackster[tsId].
empty())
355 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"trackster Id:\t" << tsId <<
"\tST id:\t-1" 361 for (
const auto& stId : sTIndices) {
362 float invSTEnergyWeight = 0.f;
364 const unsigned int STNumberOfLCs = tssInSimTrackster[stId].lcs_and_fractions.size();
365 if (STNumberOfLCs == 0)
369 int tsWithMaxEnergyInST = -1;
371 float maxEnergyTSinST = 0.f;
372 float STenergy = tssInSimTrackster[stId].energy;
374 float STEnergyFractionInTS = 0.f;
375 for (
const auto& ts : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
376 if (ts.second.first > maxEnergyTSinST) {
377 maxEnergyTSinST = ts.second.first;
378 tsWithMaxEnergyInST = ts.first;
382 STEnergyFractionInTS = maxEnergyTSinST / STenergy;
384 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
385 << std::setw(12) <<
"simTrackster\t" << std::setw(15) <<
"st total energy\t" << std::setw(15)
386 <<
"stEnergyOnLayer\t" << std::setw(14) <<
"STNhitsOnLayer\t" << std::setw(18) <<
"tsWithMaxEnergyInST\t" 387 << std::setw(15) <<
"maxEnergyTSinST\t" << std::setw(20) <<
"STEnergyFractionInTS" 389 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
390 << std::setw(12) << stId <<
"\t" << std::setw(15) << simTracksters[stId].raw_energy() <<
"\t" << std::setw(15)
391 << STenergy <<
"\t" << std::setw(14) << STNumberOfLCs <<
"\t" << std::setw(18) << tsWithMaxEnergyInST <<
"\t" 392 << std::setw(15) << maxEnergyTSinST <<
"\t" << std::setw(20) << STEnergyFractionInTS <<
"\n";
396 for (
auto const& lcaf : tssInSimTrackster[stId].lcs_and_fractions) {
399 invSTEnergyWeight = 1.f / invSTEnergyWeight;
401 for (
unsigned int i = 0;
i < STNumberOfLCs; ++
i) {
402 auto& st_lcId = tssInSimTrackster[stId].lcs_and_fractions[
i].first;
403 auto& st_lcFraction = tssInSimTrackster[stId].lcs_and_fractions[
i].second;
405 if (st_lcFraction == 0.
f)
408 for (
auto& tsPair : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
409 unsigned int tsId = tsPair.first;
411 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
412 const auto tsFraction = 1.f / tracksters[tsId].vertex_multiplicity(
i);
414 tsPair.second.second +=
415 (tsFraction - st_lcFraction) * (tsFraction - st_lcFraction) * lcEnergyWeight * invSTEnergyWeight;
417 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
418 <<
"STLCId:\t" << (uint32_t)st_lcId <<
"\tTracksterId:\t" << tsId <<
"\t" 419 <<
"tsFraction, stFraction:\t" << tsFraction <<
", " << st_lcFraction <<
"\t" 420 <<
"lcEnergyWeight:\t" << lcEnergyWeight <<
"\t" 421 <<
"current score:\t" << tsPair.second.second <<
"\t" 422 <<
"invSTEnergyWeight:\t" << invSTEnergyWeight <<
"\n";
428 if (tssInSimTrackster[stId].tracksterIdToEnergyAndScore.empty())
429 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl") <<
"ST Id:\t" << stId <<
"\tTS id:\t-1 " 432 for (
const auto& tsPair : tssInSimTrackster[stId].tracksterIdToEnergyAndScore) {
433 LogDebug(
"TSToSimTSAssociatorByEnergyScoreImpl")
434 <<
"ST Id: \t" << stId <<
"\t TS id: \t" << tsPair.first <<
"\t score \t" << tsPair.second.second
435 <<
"\t shared energy:\t" << tsPair.second.first <<
"\t shared energy fraction:\t" 436 << (tsPair.second.first / STenergy) <<
"\n";
440 return {stsInTrackster, tssInSimTrackster};
const std::unordered_map< DetId, const unsigned int > * hitMap_
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< std::vector< std::pair< unsigned int, std::pair< float, float > > > > tracksterToSimTrackster
std::vector< ticl::simTracksterOnLayer > simTracksterToTrackster
std::vector< const HGCRecHit * > hits_
Power< A, B >::type pow(const A &a, const B &b)