25 const auto& tracksters = *tCH.
product();
27 const auto& simClusters = *sCCH.
product();
28 auto nTracksters = tracksters.size();
32 auto nSimClusters = simClusters.size();
33 std::vector<size_t> sCIndices;
34 for (
unsigned int scId = 0; scId < nSimClusters; ++scId) {
36 simClusters[scId].g4Tracks()[0].eventId().bunchCrossing() != 0)) {
37 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
38 <<
"Excluding SimCluster from event: " << simClusters[scId].g4Tracks()[0].eventId().event()
39 <<
" with BX: " << simClusters[scId].g4Tracks()[0].eventId().bunchCrossing() << std::endl;
42 sCIndices.emplace_back(scId);
44 nSimClusters = sCIndices.size();
50 tssInSimCluster.resize(nSimClusters);
51 for (
unsigned int i = 0;
i < nSimClusters; ++
i) {
52 tssInSimCluster[
i].simClusterId =
i;
53 tssInSimCluster[
i].energy = 0.f;
54 tssInSimCluster[
i].hits_and_fractions.clear();
58 std::unordered_map<DetId, std::vector<ticl::detIdInfoInCluster>> detIdToSimClusterId_Map;
59 for (
const auto& scId : sCIndices) {
60 const auto& hits_and_fractions = simClusters[scId].hits_and_fractions();
61 for (
const auto& it_haf : hits_and_fractions) {
62 const auto hitid = it_haf.first;
63 const auto itcheck =
hitMap_->find(hitid);
64 if (itcheck !=
hitMap_->end()) {
65 const auto hit_find_it = detIdToSimClusterId_Map.find(hitid);
66 if (hit_find_it == detIdToSimClusterId_Map.end()) {
67 detIdToSimClusterId_Map[hitid] = std::vector<ticl::detIdInfoInCluster>();
69 detIdToSimClusterId_Map[hitid].emplace_back(scId, it_haf.second);
72 tssInSimCluster[scId].energy += it_haf.second *
hit->energy();
73 tssInSimCluster[scId].hits_and_fractions.emplace_back(hitid, it_haf.second);
79 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
80 <<
"tssInSimCluster INFO (Only SimCluster filled at the moment)" << std::endl;
81 for (
size_t sc = 0; sc < tssInSimCluster.size(); ++sc) {
82 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"For SimCluster Idx: " << sc <<
" we have: " << std::endl;
83 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
84 <<
"\tSimClusterIdx:\t" << tssInSimCluster[sc].simClusterId << std::endl;
85 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\tEnergy:\t" << tssInSimCluster[sc].energy << std::endl;
86 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\t# of clusters:\t" <<
layerClusters.size() << std::endl;
87 double tot_energy = 0.;
88 for (
auto const& haf : tssInSimCluster[sc].hits_and_fractions) {
90 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\tHits/fraction/energy: " << (uint32_t)haf.first <<
"/" 91 << haf.second <<
"/" << haf.second *
hit->energy() << std::endl;
92 tot_energy += haf.second *
hit->energy();
94 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\tTot Sum haf: " << tot_energy << std::endl;
95 for (
auto const& ts : tssInSimCluster[sc].tracksterIdToEnergyAndScore) {
96 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
97 <<
"\ttsIdx/energy/score: " << ts.first <<
"/" << ts.second.first <<
"/" << ts.second.second << std::endl;
101 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"detIdToSimClusterId_Map INFO" << std::endl;
102 for (
auto const&
detId : detIdToSimClusterId_Map) {
104 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
105 <<
"For detId: " << (uint32_t)
detId.first
106 <<
" we have found the following connections with SimClusters:" << std::endl;
107 for (
auto const& sc :
detId.second) {
108 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
109 <<
"\tSimCluster Id: " << sc.clusterId <<
" with fraction: " << sc.fraction
110 <<
" and energy: " << sc.fraction *
hit->energy() << std::endl;
116 std::unordered_map<DetId, std::vector<ticl::detIdInfoInCluster>> detIdToLayerClusterId_Map;
121 scsInTrackster.resize(nTracksters);
123 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
124 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
125 const auto lcId = tracksters[tsId].vertices(
i);
126 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
128 const std::vector<std::pair<DetId, float>>& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
129 unsigned int numberOfHitsInLC = hits_and_fractions.size();
131 for (
unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
132 const auto rh_detid = hits_and_fractions[hitId].first;
133 const auto rhFraction = hits_and_fractions[hitId].second;
135 const auto hit_find_in_LC = detIdToLayerClusterId_Map.find(rh_detid);
136 if (hit_find_in_LC == detIdToLayerClusterId_Map.end()) {
137 detIdToLayerClusterId_Map[rh_detid] = std::vector<ticl::detIdInfoInCluster>();
139 detIdToLayerClusterId_Map[rh_detid].emplace_back(lcId, rhFraction);
141 const auto hit_find_in_SC = detIdToSimClusterId_Map.find(rh_detid);
143 if (hit_find_in_SC != detIdToSimClusterId_Map.end()) {
144 const auto itcheck =
hitMap_->find(rh_detid);
149 for (
const auto&
h : hit_find_in_SC->second) {
152 tssInSimCluster[
h.clusterId].tracksterIdToEnergyAndScore[tsId].first +=
153 lcFractionInTs *
h.fraction *
hit->energy();
155 scsInTrackster[tsId].emplace_back(
h.clusterId, 0.f);
163 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
164 for (
const auto& lcId : tracksters[tsId].
vertices()) {
165 const auto& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
166 unsigned int numberOfHitsInLC = hits_and_fractions.size();
176 std::vector<int> hitsToSimClusterId(numberOfHitsInLC);
178 int maxSCId_byNumberOfHits = -1;
180 unsigned int maxSCNumberOfHitsInLC = 0;
182 int maxSCId_byEnergy = -1;
184 float maxEnergySharedLCandSC = 0.f;
186 float energyFractionOfLCinSC = 0.f;
188 float energyFractionOfSCinLC = 0.f;
189 std::unordered_map<unsigned, unsigned> occurrencesSCinLC;
190 unsigned int numberOfNoiseHitsInLC = 0;
191 std::unordered_map<unsigned, float> SCEnergyInLC;
193 for (
unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
194 const auto rh_detid = hits_and_fractions[hitId].first;
195 const auto rhFraction = hits_and_fractions[hitId].second;
197 const auto hit_find_in_SC = detIdToSimClusterId_Map.find(rh_detid);
205 if (rhFraction == 0.) {
206 hitsToSimClusterId[hitId] = -2;
209 if (hit_find_in_SC == detIdToSimClusterId_Map.end()) {
210 hitsToSimClusterId[hitId] -= 1;
212 const auto itcheck =
hitMap_->find(rh_detid);
214 auto maxSCEnergyInLC = 0.f;
217 for (
const auto&
h : hit_find_in_SC->second) {
218 SCEnergyInLC[
h.clusterId] +=
h.fraction *
hit->energy();
221 if (SCEnergyInLC[
h.clusterId] > maxSCEnergyInLC) {
222 maxSCEnergyInLC = SCEnergyInLC[
h.clusterId];
223 maxSCId =
h.clusterId;
226 hitsToSimClusterId[hitId] = maxSCId;
230 for (
const auto&
c : hitsToSimClusterId) {
232 numberOfNoiseHitsInLC++;
234 occurrencesSCinLC[
c]++;
238 for (
const auto&
c : occurrencesSCinLC) {
239 if (
c.second > maxSCNumberOfHitsInLC) {
240 maxSCId_byNumberOfHits =
c.first;
241 maxSCNumberOfHitsInLC =
c.second;
245 for (
const auto&
c : SCEnergyInLC) {
246 if (
c.second > maxEnergySharedLCandSC) {
247 maxSCId_byEnergy =
c.first;
248 maxEnergySharedLCandSC =
c.second;
252 float totalSCEnergyOnLayer = 0.f;
253 if (maxSCId_byEnergy >= 0) {
254 totalSCEnergyOnLayer = tssInSimCluster[maxSCId_byEnergy].energy;
255 energyFractionOfSCinLC = maxEnergySharedLCandSC / totalSCEnergyOnLayer;
256 if (tracksters[tsId].raw_energy() > 0.
f) {
257 energyFractionOfLCinSC = maxEnergySharedLCandSC / tracksters[tsId].raw_energy();
261 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
262 << std::setw(12) <<
"TracksterID:\t" << std::setw(12) <<
"layerCluster\t" << std::setw(10) <<
"lc energy\t" 263 << std::setw(5) <<
"nhits\t" << std::setw(12) <<
"noise hits\t" << std::setw(22) <<
"maxSCId_byNumberOfHits\t" 264 << std::setw(8) <<
"nhitsSC\t" << std::setw(13) <<
"maxSCId_byEnergy\t" << std::setw(20)
265 <<
"maxEnergySharedLCandSC\t" << std::setw(22) <<
"totalSCEnergyOnLayer\t" << std::setw(22)
266 <<
"energyFractionOfLCinSC\t" << std::setw(25) <<
"energyFractionOfSCinLC\t" 268 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
269 << std::setw(12) << tsId <<
"\t" << std::setw(12) << lcId <<
"\t" << std::setw(10)
270 << tracksters[tsId].raw_energy() <<
"\t" << std::setw(5) << numberOfHitsInLC <<
"\t" << std::setw(12)
271 << numberOfNoiseHitsInLC <<
"\t" << std::setw(22) << maxSCId_byNumberOfHits <<
"\t" << std::setw(8)
272 << maxSCNumberOfHitsInLC <<
"\t" << std::setw(13) << maxSCId_byEnergy <<
"\t" << std::setw(20)
273 << maxEnergySharedLCandSC <<
"\t" << std::setw(22) << totalSCEnergyOnLayer <<
"\t" << std::setw(22)
274 << energyFractionOfLCinSC <<
"\t" << std::setw(25) << energyFractionOfSCinLC <<
"\n";
278 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
279 <<
"Improved tssInSimCluster INFO (Now containing the linked tracksters id and energy - score still empty)" 281 for (
size_t sc = 0; sc < tssInSimCluster.size(); ++sc) {
282 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"For SimCluster Idx: " << sc <<
" we have: " << std::endl;
283 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
284 <<
" SimClusterIdx: " << tssInSimCluster[sc].simClusterId << std::endl;
285 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\tEnergy:\t" << tssInSimCluster[sc].energy << std::endl;
286 double tot_energy = 0.;
287 for (
auto const& haf : tssInSimCluster[sc].hits_and_fractions) {
289 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\tHits/fraction/energy: " << (uint32_t)haf.first <<
"/" 290 << haf.second <<
"/" << haf.second *
hit->energy() << std::endl;
291 tot_energy += haf.second *
hit->energy();
293 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"\tTot Sum haf: " << tot_energy << std::endl;
294 for (
auto const& ts : tssInSimCluster[sc].tracksterIdToEnergyAndScore) {
295 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
296 <<
"\ttsIdx/energy/score: " << ts.first <<
"/" << ts.second.first <<
"/" << ts.second.second << std::endl;
300 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"Improved detIdToSimClusterId_Map INFO" << std::endl;
301 for (
auto const& sc : detIdToSimClusterId_Map) {
303 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
304 <<
"For detId: " << (uint32_t)sc.first
305 <<
" we have found the following connections with SimClusters:" << std::endl;
306 for (
auto const& sclu : sc.second) {
307 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
308 <<
" SimCluster Id: " << sclu.clusterId <<
" with fraction: " << sclu.fraction
309 <<
" and energy: " << sclu.fraction *
hit->energy() << std::endl;
316 for (
unsigned int tsId = 0; tsId < nTracksters; ++tsId) {
321 scsInTrackster[tsId].erase(
last, scsInTrackster[tsId].
end());
325 if (tracksters[tsId].raw_energy() == 0. && !scsInTrackster[tsId].
empty()) {
326 for (
auto& scPair : scsInTrackster[tsId]) {
328 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
329 <<
"TracksterId:\t " << tsId <<
"\tSC id:\t" << scPair.first <<
"\tscore\t " << scPair.second <<
"\n";
334 float invTracksterEnergyWeight = 0.f;
335 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
336 const auto lcId = tracksters[tsId].vertices(
i);
337 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
339 const auto& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
341 for (
auto const& haf : hits_and_fractions) {
343 invTracksterEnergyWeight +=
344 (lcFractionInTs * haf.second *
hit->energy()) * (lcFractionInTs * haf.second *
hit->energy());
347 invTracksterEnergyWeight = 1.f / invTracksterEnergyWeight;
349 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
350 const auto lcId = tracksters[tsId].vertices(
i);
351 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
353 const auto& hits_and_fractions =
layerClusters[lcId].hitsAndFractions();
354 unsigned int numberOfHitsInLC = hits_and_fractions.size();
355 for (
unsigned int i = 0;
i < numberOfHitsInLC; ++
i) {
356 DetId rh_detid = hits_and_fractions[
i].first;
357 float rhFraction = hits_and_fractions[
i].second * lcFractionInTs;
359 const bool hitWithSC = (detIdToSimClusterId_Map.find(rh_detid) != detIdToSimClusterId_Map.end());
361 const auto itcheck =
hitMap_->find(rh_detid);
363 float hitEnergyWeight =
hit->energy() *
hit->energy();
365 for (
auto& scPair : scsInTrackster[tsId]) {
366 float scFraction = 0.f;
368 const auto findHitIt =
std::find(detIdToSimClusterId_Map[rh_detid].
begin(),
369 detIdToSimClusterId_Map[rh_detid].
end(),
371 if (findHitIt != detIdToSimClusterId_Map[rh_detid].
end())
372 scFraction = findHitIt->fraction;
375 (rhFraction - scFraction) * (rhFraction - scFraction) * hitEnergyWeight * invTracksterEnergyWeight;
377 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
378 <<
"rh_detid:\t" << (uint32_t)rh_detid <<
"\ttracksterId:\t" << tsId <<
"\t" 379 <<
"rhfraction,scFraction:\t" << rhFraction <<
", " << scFraction <<
"\t" 380 <<
"hitEnergyWeight:\t" << hitEnergyWeight <<
"\t" 381 <<
"current score:\t" << scPair.second <<
"\t" 382 <<
"invTracksterEnergyWeight:\t" << invTracksterEnergyWeight <<
"\n";
389 if (scsInTrackster[tsId].
empty())
390 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"trackster Id:\t" << tsId <<
"\tSC id:\t-1" 396 for (
const auto& scId : sCIndices) {
397 float invSCEnergyWeight = 0.f;
399 const unsigned int SCNumberOfHits = tssInSimCluster[scId].hits_and_fractions.size();
400 if (SCNumberOfHits == 0)
403 int tsWithMaxEnergyInSC = -1;
405 float maxEnergyTSinSC = 0.f;
406 float SCenergy = tssInSimCluster[scId].energy;
408 float SCEnergyFractionInTS = 0.f;
409 for (
const auto& ts : tssInSimCluster[scId].tracksterIdToEnergyAndScore) {
410 if (ts.second.first > maxEnergyTSinSC) {
411 maxEnergyTSinSC = ts.second.first;
412 tsWithMaxEnergyInSC = ts.first;
416 SCEnergyFractionInTS = maxEnergyTSinSC / SCenergy;
418 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
419 << std::setw(12) <<
"simcluster\t" << std::setw(15) <<
"sc total energy\t" << std::setw(15)
420 <<
"scEnergyOnLayer\t" << std::setw(14) <<
"SCNhitsOnLayer\t" << std::setw(18) <<
"tsWithMaxEnergyInSC\t" 421 << std::setw(15) <<
"maxEnergyTSinSC\t" << std::setw(20) <<
"SCEnergyFractionInTS" 423 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
424 << std::setw(12) << scId <<
"\t" << std::setw(15) << simClusters[scId].energy() <<
"\t" << std::setw(15)
425 << SCenergy <<
"\t" << std::setw(14) << SCNumberOfHits <<
"\t" << std::setw(18) << tsWithMaxEnergyInSC <<
"\t" 426 << std::setw(15) << maxEnergyTSinSC <<
"\t" << std::setw(20) << SCEnergyFractionInTS <<
"\n";
429 for (
auto const& haf : tssInSimCluster[scId].hits_and_fractions) {
431 invSCEnergyWeight +=
std::pow(haf.second *
hit->energy(), 2);
433 invSCEnergyWeight = 1.f / invSCEnergyWeight;
435 for (
unsigned int i = 0;
i < SCNumberOfHits; ++
i) {
436 auto& sc_hitDetId = tssInSimCluster[scId].hits_and_fractions[
i].first;
437 auto& scFraction = tssInSimCluster[scId].hits_and_fractions[
i].second;
439 bool hitWithLC =
false;
440 if (scFraction == 0.
f)
442 const auto hit_find_in_LC = detIdToLayerClusterId_Map.find(sc_hitDetId);
443 if (hit_find_in_LC != detIdToLayerClusterId_Map.end())
445 const auto itcheck =
hitMap_->find(sc_hitDetId);
447 float hitEnergyWeight =
hit->energy() *
hit->energy();
448 for (
auto& tsPair : tssInSimCluster[scId].tracksterIdToEnergyAndScore) {
449 unsigned int tsId = tsPair.first;
450 float tsFraction = 0.f;
452 for (
unsigned int i = 0;
i < tracksters[tsId].vertices().size(); ++
i) {
453 const auto lcId = tracksters[tsId].vertices(
i);
454 const auto lcFractionInTs = 1.f / tracksters[tsId].vertex_multiplicity(
i);
457 const auto findHitIt =
std::find(detIdToLayerClusterId_Map[sc_hitDetId].
begin(),
458 detIdToLayerClusterId_Map[sc_hitDetId].
end(),
460 if (findHitIt != detIdToLayerClusterId_Map[sc_hitDetId].
end())
461 tsFraction = findHitIt->fraction * lcFractionInTs;
463 tsPair.second.second +=
464 (tsFraction - scFraction) * (tsFraction - scFraction) * hitEnergyWeight * invSCEnergyWeight;
466 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
467 <<
"SCDetId:\t" << (uint32_t)sc_hitDetId <<
"\tTracksterId:\t" << tsId <<
"\t" 468 <<
"tsFraction, scFraction:\t" << tsFraction <<
", " << scFraction <<
"\t" 469 <<
"hitEnergyWeight:\t" << hitEnergyWeight <<
"\t" 470 <<
"current score:\t" << tsPair.second.second <<
"\t" 471 <<
"invSCEnergyWeight:\t" << invSCEnergyWeight <<
"\n";
477 if (tssInSimCluster[scId].tracksterIdToEnergyAndScore.empty())
478 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl") <<
"SC Id:\t" << scId <<
"\tTS id:\t-1 " 481 for (
const auto& tsPair : tssInSimCluster[scId].tracksterIdToEnergyAndScore) {
482 LogDebug(
"TSToSCAssociatorByEnergyScoreImpl")
483 <<
"SC Id: \t" << scId <<
"\t TS id: \t" << tsPair.first <<
"\t score \t" << tsPair.second.second
484 <<
"\t shared energy:\t" << tsPair.second.first <<
"\t shared energy fraction:\t" 485 << (tsPair.second.first / SCenergy) <<
"\n";
489 return {scsInTrackster, tssInSimCluster};
for(int i=first, nt=offsets[nh];i< nt;i+=gridDim.x *blockDim.x)
std::vector< const HGCRecHit * > hits_
T const * product() const
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
const bool hardScatterOnly_
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
def unique(seq, keepstr=True)
const std::unordered_map< DetId, const unsigned int > * hitMap_
std::vector< std::vector< std::pair< unsigned int, float > > > tracksterToSimCluster
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Power< A, B >::type pow(const A &a, const B &b)
std::vector< ticl::simClusterOnBLayer > simClusterToTrackster