CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
TracksterToSimTracksterAssociatorProducer.cc
Go to the documentation of this file.
1 // Author: Felice Pantaleo, felice.pantaleo@cern.ch 06/2024
2 
5 
7  : recoTracksterCollectionToken_(
8  consumes<std::vector<ticl::Trackster>>(pset.getParameter<edm::InputTag>("tracksters"))),
9  simTracksterCollectionToken_(
10  consumes<std::vector<ticl::Trackster>>(pset.getParameter<edm::InputTag>("simTracksters"))),
11  layerClustersCollectionToken_(
12  consumes<std::vector<reco::CaloCluster>>(pset.getParameter<edm::InputTag>("layerClusters"))),
13  LayerClusterToTracksterMapToken_(
14  consumes<
15  ticl::AssociationMap<ticl::mapWithFraction, std::vector<reco::CaloCluster>, std::vector<ticl::Trackster>>>(
16  pset.getParameter<edm::InputTag>("tracksterMap"))),
17  LayerClusterToSimTracksterMapToken_(
18  consumes<
19  ticl::AssociationMap<ticl::mapWithFraction, std::vector<reco::CaloCluster>, std::vector<ticl::Trackster>>>(
20  pset.getParameter<edm::InputTag>("simTracksterMap"))) {
21  produces<
23  "tracksterToSimTracksterMap");
24  produces<
26  "simTracksterToTracksterMap");
27 }
28 
30 
33  const edm::EventSetup& iSetup) const {
34  edm::Handle<std::vector<ticl::Trackster>> recoTrackstersHandle;
35  iEvent.getByToken(recoTracksterCollectionToken_, recoTrackstersHandle);
36  const auto& recoTracksters = *recoTrackstersHandle;
37 
38  edm::Handle<std::vector<ticl::Trackster>> simTrackstersHandle;
39  iEvent.getByToken(simTracksterCollectionToken_, simTrackstersHandle);
40  const auto& simTracksters = *simTrackstersHandle;
41 
42  edm::Handle<std::vector<reco::CaloCluster>> layerClustersHandle;
43  iEvent.getByToken(layerClustersCollectionToken_, layerClustersHandle);
44  const auto& layerClusters = *layerClustersHandle;
45 
47  layerClusterToTracksterMapHandle;
48  iEvent.getByToken(LayerClusterToTracksterMapToken_, layerClusterToTracksterMapHandle);
49  const auto& layerClusterToTracksterMap = *layerClusterToTracksterMapHandle;
50 
52  layerClusterToSimTracksterMapHandle;
53  iEvent.getByToken(LayerClusterToSimTracksterMapToken_, layerClusterToSimTracksterMapHandle);
54  const auto& layerClusterToSimTracksterMap = *layerClusterToSimTracksterMapHandle;
55 
56  auto tracksterToSimTracksterMap = std::make_unique<
58  recoTrackstersHandle, simTrackstersHandle, iEvent);
59  auto simTracksterToTracksterMap = std::make_unique<
61  simTrackstersHandle, recoTrackstersHandle, iEvent);
62 
63  for (unsigned int tracksterIndex = 0; tracksterIndex < recoTracksters.size(); ++tracksterIndex) {
64  const auto& recoTrackster = recoTracksters[tracksterIndex];
65  edm::Ref<std::vector<ticl::Trackster>> recoTracksterRef(recoTrackstersHandle, tracksterIndex);
66  const auto& layerClustersIds = recoTrackster.vertices();
67  float recoToSimScoresDenominator = 0.f;
68  ticl::mapWithFraction layerClusterToAssociatedSimTracksterMap(layerClustersIds.size());
69  std::vector<unsigned int> associatedSimTracksterIndices;
70  for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
71  unsigned int layerClusterId = layerClustersIds[i];
72  const auto& layerCluster = layerClusters[layerClusterId];
73  float recoFraction = 1.f / recoTrackster.vertex_multiplicity(i);
74  float squaredRecoFraction = recoFraction * recoFraction;
75  float squaredLayerClusterEnergy = layerCluster.energy() * layerCluster.energy();
76  recoToSimScoresDenominator += squaredLayerClusterEnergy * squaredRecoFraction;
77  const auto& simTracksterVec = layerClusterToSimTracksterMap[layerClusterId];
78  for (const auto& [simTracksterIndex, simSharedEnergy] : simTracksterVec) {
79  layerClusterToAssociatedSimTracksterMap[i].push_back({simTracksterIndex, simSharedEnergy});
80  associatedSimTracksterIndices.push_back(simTracksterIndex);
81  }
82  }
83 
84  // Keep only unique associatedSimTracksterIndices
85  std::sort(associatedSimTracksterIndices.begin(), associatedSimTracksterIndices.end());
86  associatedSimTracksterIndices.erase(
87  std::unique(associatedSimTracksterIndices.begin(), associatedSimTracksterIndices.end()),
88  associatedSimTracksterIndices.end());
89 
90  // Add missing sim tracksters with 0 shared energy to layerClusterToAssociatedSimTracksterMap
91  for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
92  unsigned int layerClusterId = layerClustersIds[i];
93  const auto& simTracksterVec = layerClusterToSimTracksterMap[layerClusterId];
94  for (unsigned int simTracksterIndex : associatedSimTracksterIndices) {
95  if (std::find_if(simTracksterVec.begin(), simTracksterVec.end(), [simTracksterIndex](const auto& pair) {
96  return pair.first == simTracksterIndex;
97  }) == simTracksterVec.end()) {
98  layerClusterToAssociatedSimTracksterMap[i].push_back({simTracksterIndex, 0.f});
99  }
100  }
101  }
102 
103  const float invDenominator = 1.f / recoToSimScoresDenominator;
104 
105  for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
106  unsigned int layerClusterId = layerClustersIds[i];
107  const auto& layerCluster = layerClusters[layerClusterId];
108  float recoFraction = 1.f / recoTrackster.vertex_multiplicity(i);
109  float squaredRecoFraction = recoFraction * recoFraction;
110  float squaredLayerClusterEnergy = layerCluster.energy() * layerCluster.energy();
111  float recoSharedEnergy = layerCluster.energy() * recoFraction;
112  float invLayerClusterEnergy = 1.f / layerCluster.energy();
113  const auto& simTracksterVec = layerClusterToAssociatedSimTracksterMap[i];
114  for (const auto& [simTracksterIndex, simSharedEnergy] : simTracksterVec) {
115  edm::Ref<std::vector<ticl::Trackster>> simTracksterRef(simTrackstersHandle, simTracksterIndex);
116  float sharedEnergy = std::min(simSharedEnergy, recoSharedEnergy);
117  float simFraction = simSharedEnergy * invLayerClusterEnergy;
118  float score = invDenominator *
119  std::min(squaredRecoFraction, (recoFraction - simFraction) * (recoFraction - simFraction)) *
120  squaredLayerClusterEnergy;
121  tracksterToSimTracksterMap->insert(recoTracksterRef, simTracksterRef, sharedEnergy, score);
122  }
123  }
124  }
125 
126  for (unsigned int tracksterIndex = 0; tracksterIndex < simTracksters.size(); ++tracksterIndex) {
127  const auto& simTrackster = simTracksters[tracksterIndex];
128  edm::Ref<std::vector<ticl::Trackster>> simTracksterRef(simTrackstersHandle, tracksterIndex);
129  const auto& layerClustersIds = simTrackster.vertices();
130  float simToRecoScoresDenominator = 0.f;
131  ticl::mapWithFraction layerClusterToAssociatedTracksterMap(layerClustersIds.size());
132  std::vector<unsigned int> associatedRecoTracksterIndices;
133  for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
134  unsigned int layerClusterId = layerClustersIds[i];
135  const auto& layerCluster = layerClusters[layerClusterId];
136  float simFraction = 1.f / simTrackster.vertex_multiplicity(i);
137  float squaredSimFraction = simFraction * simFraction;
138  float squaredLayerClusterEnergy = layerCluster.energy() * layerCluster.energy();
139  simToRecoScoresDenominator += squaredLayerClusterEnergy * squaredSimFraction;
140  const auto& recoTracksterVec = layerClusterToTracksterMap[layerClusterId];
141  for (const auto& [recoTracksterIndex, recoSharedEnergy] : recoTracksterVec) {
142  layerClusterToAssociatedTracksterMap[i].push_back({recoTracksterIndex, recoSharedEnergy});
143  associatedRecoTracksterIndices.push_back(recoTracksterIndex);
144  }
145  }
146  // keep only unique associatedRecoTracksterIndices
147  std::sort(associatedRecoTracksterIndices.begin(), associatedRecoTracksterIndices.end());
148  associatedRecoTracksterIndices.erase(
149  std::unique(associatedRecoTracksterIndices.begin(), associatedRecoTracksterIndices.end()),
150  associatedRecoTracksterIndices.end());
151  // for each layer cluster, loop over associatedRecoTracksterIndices and add the missing reco tracksters with 0 shared energy to layerClusterToAssociatedTracksterMap
152  for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
153  unsigned int layerClusterId = layerClustersIds[i];
154  const auto& recoTracksterVec = layerClusterToTracksterMap[layerClusterId];
155  for (unsigned int recoTracksterIndex : associatedRecoTracksterIndices) {
156  if (std::find_if(recoTracksterVec.begin(), recoTracksterVec.end(), [recoTracksterIndex](const auto& pair) {
157  return pair.first == recoTracksterIndex;
158  }) == recoTracksterVec.end()) {
159  layerClusterToAssociatedTracksterMap[i].push_back({recoTracksterIndex, 0.f});
160  }
161  }
162  }
163 
164  const float invDenominator = 1.f / simToRecoScoresDenominator;
165 
166  for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
167  unsigned int layerClusterId = layerClustersIds[i];
168  const auto& layerCluster = layerClusters[layerClusterId];
169  float simFraction = 1.f / simTrackster.vertex_multiplicity(i);
170  float squaredSimFraction = simFraction * simFraction;
171  float squaredLayerClusterEnergy = layerCluster.energy() * layerCluster.energy();
172  float simSharedEnergy = layerCluster.energy() * simFraction;
173  float invLayerClusterEnergy = 1.f / layerCluster.energy();
174  const auto& recoTracksterVec = layerClusterToAssociatedTracksterMap[i];
175  for (const auto& [recoTracksterIndex, recoSharedEnergy] : recoTracksterVec) {
176  edm::Ref<std::vector<ticl::Trackster>> recoTracksterRef(recoTrackstersHandle, recoTracksterIndex);
177  float sharedEnergy = std::min(recoSharedEnergy, simSharedEnergy);
178  float recoFraction = recoSharedEnergy * invLayerClusterEnergy;
179  float score = invDenominator *
180  std::min(squaredSimFraction, (simFraction - recoFraction) * (simFraction - recoFraction)) *
181  squaredLayerClusterEnergy;
182  simTracksterToTracksterMap->insert(tracksterIndex, recoTracksterIndex, sharedEnergy, score);
183  }
184  }
185  }
186  tracksterToSimTracksterMap->sort(true);
187  simTracksterToTracksterMap->sort(true);
188  iEvent.put(std::move(tracksterToSimTracksterMap), "tracksterToSimTracksterMap");
189  iEvent.put(std::move(simTracksterToTracksterMap), "simTracksterToTracksterMap");
190 }
191 
194  desc.add<edm::InputTag>("tracksters", edm::InputTag("trackstersProducer"));
195  desc.add<edm::InputTag>("simTracksters", edm::InputTag("simTrackstersProducer"));
196  desc.add<edm::InputTag>("layerClusters", edm::InputTag("hgcalMergeLayerClusters"));
197  desc.add<edm::InputTag>("tracksterMap", edm::InputTag("tracksterAssociatorProducer"));
198  desc.add<edm::InputTag>("simTracksterMap", edm::InputTag("simTracksterAssociatorProducer"));
199  descriptions.add("tracksterToSimTracksterAssociatorProducer", desc);
200 }
201 
edm::EDGetTokenT< std::vector< reco::CaloCluster > > layerClustersCollectionToken_
edm::EDGetTokenT< std::vector< ticl::Trackster > > simTracksterCollectionToken_
std::vector< std::vector< std::pair< unsigned int, float > >> mapWithFraction
int iEvent
Definition: GenABIO.cc:224
def unique(seq, keepstr=True)
Definition: tier0.py:24
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
edm::EDGetTokenT< ticl::AssociationMap< ticl::mapWithFraction, std::vector< reco::CaloCluster >, std::vector< ticl::Trackster > > > LayerClusterToTracksterMapToken_
edm::EDGetTokenT< ticl::AssociationMap< ticl::mapWithFraction, std::vector< reco::CaloCluster >, std::vector< ticl::Trackster > > > LayerClusterToSimTracksterMapToken_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override
void add(std::string const &label, ParameterSetDescription const &psetDescription)
float sharedEnergy(reco::CaloCluster const &clu1, reco::CaloCluster const &clu2, EcalRecHitCollection const &barrelRecHits, EcalRecHitCollection const &endcapRecHits)
fixed size matrix
HLT enums.
Definition: Common.h:10
edm::EDGetTokenT< std::vector< ticl::Trackster > > recoTracksterCollectionToken_
def move(src, dest)
Definition: eostools.py:511