Process the event and put new data products into it.
132 {
133 if (!event.
exists(input_track_coll_name_, input_tracks_passname_)) {
134 ldmx_log(error) << "Unable to find (one) collection named "
135 << input_track_coll_name_ << "_" << input_tracks_passname_;
136 return;
137 }
138 if (!event.
exists(input_ecal_coll_name_, input_ecal_passname_)) {
139 ldmx_log(error) << "Unable to find (one) collection named "
140 << input_ecal_coll_name_ << "_" << input_ecal_passname_;
141 return;
142 }
143 if (!event.
exists(input_hcal_coll_name_, input_hcal_passname_)) {
144 ldmx_log(error) << "Unable to find (one) collection named "
145 << input_hcal_coll_name_ << "_" << input_hcal_passname_;
146 return;
147 }
148
150 input_hcal_coll_name_, input_hcal_passname_);
152 input_track_coll_name_, input_tracks_passname_);
153
154 const auto ecal_clusters =
155 use_existing_ecal_clusters_
156 ? getEcalClusters(event, input_ecal_coll_name_, input_ecal_passname_)
157 : event.getCollection<ldmx::CaloCluster>(input_ecal_coll_name_,
158 input_ecal_passname_);
159
160 std::vector<ldmx::PFCandidate> pf_cands;
161
162 if (!single_particle_) {
163
164
165
166
167
168
169
170
171
172
173
174
175
176 std::map<int, std::vector<int> > tk_calo_map;
177 std::map<int, std::vector<int> > calo_tk_map;
178 std::map<std::pair<int, int>, float> tk_em_dist;
179
180 for (int i = 0; i < tracks.size(); i++) {
181 const auto& tk = tracks[i];
182 const std::vector<float> xyz = tk.getPosition();
183 const std::vector<double> pxyz = tk.getMomentum();
184 const float p = sqrt(pow(pxyz[0], 2) + pow(pxyz[1], 2) + pow(pxyz[2], 2));
185
186 for (int j = 0; j < ecal_clusters.size(); j++) {
187 const auto& ecal = ecal_clusters[j];
188
189 const float ecal_clus_z = ecal.getCentroidZ();
190 const float tk_x_at_clus =
191 xyz[0] +
192 pxyz[0] / pxyz[2] * (ecal_clus_z - xyz[2]);
193 const float tk_y_at_clus =
194 xyz[1] + pxyz[1] / pxyz[2] * (ecal_clus_z - xyz[2]);
195 float dist = hypot((tk_x_at_clus - ecal.getCentroidX()) /
196 std::max(1.0, ecal.getRMSX()),
197 (tk_y_at_clus - ecal.getCentroidY()) /
198 std::max(1.0, ecal.getRMSY()));
199 tk_em_dist[{i, j}] = dist;
200 bool is_match =
201 (dist < 2) && (ecal.getEnergy() > 0.3 * p &&
202 ecal.getEnergy() < 2 * p);
203
204 if (is_match) {
205 if (tk_calo_map.count(i))
206 tk_calo_map[i].push_back(j);
207 else
208 tk_calo_map[i] = {j};
209 if (calo_tk_map.count(j))
210 calo_tk_map[j].push_back(i);
211 else
212 calo_tk_map[j] = {i};
213 }
214 }
215 }
216
217
218 std::map<int, std::vector<int> > em_had_calo_map;
219 std::map<std::pair<int, int>, float> em_had_dist;
220 for (int i = 0; i < ecal_clusters.size(); i++) {
221 const auto& ecal = ecal_clusters[i];
222 for (int j = 0; j < hcal_clusters.size(); j++) {
223 const auto& hcal = hcal_clusters[j];
224
225 const float x_at_h_clus =
226 ecal.getCentroidX() +
227 ecal.getDXDZ() * (hcal.getCentroidZ() -
228 ecal.getCentroidZ());
229 const float y_at_h_clus =
230 ecal.getCentroidY() +
231 ecal.getDYDZ() * (hcal.getCentroidZ() - ecal.getCentroidZ());
232 float dist = sqrt(
233 pow(x_at_h_clus - hcal.getCentroidX(), 2) /
234 std::max(1.0, pow(hcal.getRMSX(), 2) + pow(ecal.getRMSX(), 2)) +
235 pow(y_at_h_clus - hcal.getCentroidY(), 2) /
236 std::max(1.0, pow(hcal.getRMSY(), 2) + pow(ecal.getRMSY(), 2)));
237 em_had_dist[{i, j}] = dist;
238 bool is_match = (dist < 5);
239 if (is_match) {
240 if (em_had_calo_map.count(i))
241 em_had_calo_map[i].push_back(j);
242 else
243 em_had_calo_map[i] = {j};
244 }
245 }
246 }
247
248
249
250 std::map<int, std::vector<int> > tk_had_calo_map;
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267 std::vector<bool> tk_is_em_linked(tracks.size(), false);
268 std::vector<bool> em_is_tk_linked(ecal_clusters.size(), false);
269 std::map<int, int> tk_em_pairs{};
270 for (int i = 0; i < tracks.size(); i++) {
271 if (tk_calo_map.count(i)) {
272
273 for (int em_idx : tk_calo_map[i]) {
274 if (!em_is_tk_linked[em_idx]) {
275 em_is_tk_linked[em_idx] = true;
276 tk_is_em_linked[i] = true;
277 tk_em_pairs[i] = em_idx;
278 break;
279 }
280 }
281 }
282 }
283
284
285 std::vector<bool> em_is_had_linked(ecal_clusters.size(), false);
286 std::vector<bool> had_is_em_linked(hcal_clusters.size(), false);
287 std::map<int, int> em_had_pairs{};
288 for (int i = 0; i < ecal_clusters.size(); i++) {
289 if (em_had_calo_map.count(i)) {
290
291 for (int had_idx : em_had_calo_map[i]) {
292 if (!had_is_em_linked[had_idx]) {
293 had_is_em_linked[had_idx] = true;
294 em_is_had_linked[i] = true;
295 em_had_pairs[i] = had_idx;
296 break;
297 }
298 }
299 }
300 }
301
302
303
304
305
306
307
308
309
310
311
312 for (int i = 0; i < tracks.size(); i++) {
314 fillCandTrack(cand, tracks[i]);
315
316 if (!tk_is_em_linked[i]) {
317
318 } else {
319 fillCandEMCalo(cand, ecal_clusters[tk_em_pairs[i]]);
320 if (em_is_had_linked[tk_em_pairs[i]]) {
321
322 fillCandHadCalo(cand, hcal_clusters[em_had_pairs[tk_em_pairs[i]]]);
323 }
324
325 }
326 pf_cands.push_back(cand);
327 }
328
329
330
331 for (int i = 0; i < ecal_clusters.size(); i++) {
332
333 if (em_is_tk_linked[i]) continue;
334
336 fillCandEMCalo(cand, ecal_clusters[i]);
337 if (em_is_had_linked[tk_em_pairs[i]]) {
338 fillCandHadCalo(cand, hcal_clusters[em_had_pairs[i]]);
339
340 } else {
341
342 }
343 pf_cands.push_back(cand);
344 }
345 std::vector<ldmx::PFCandidate> had_only;
346 for (int i = 0; i < hcal_clusters.size(); i++) {
347 if (had_is_em_linked[i]) continue;
349 fillCandHadCalo(cand, hcal_clusters[i]);
350
351 pf_cands.push_back(cand);
352 }
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424 } else {
425
427 int pid = 0;
428 if (tracks.size()) {
429 fillCandTrack(pf, tracks[0]);
430 pid += 1;
431 }
432 if (ecal_clusters.size()) {
433 fillCandEMCalo(pf, ecal_clusters[0]);
434 pid += 2;
435 }
436 if (hcal_clusters.size()) {
437 fillCandHadCalo(pf, hcal_clusters[0]);
438 pid += 4;
439 }
440 pf.setPID(pid);
441 pf.setEnergy(pf.getEcalEnergy() + pf.getHcalEnergy());
442 pf_cands.push_back(pf);
443 }
444
445 event.
add(output_coll_name_, pf_cands);
446}
bool exists(const std::string &name, const std::string &passName, bool unique=true) const
Check for the existence of an object or collection with the given name and pass name in the event.
void add(const std::string &collectionName, T &obj)
Adds an object to the event bus.
Represents a reconstructed particle.
Represents a simulated tracker hit in the simulation.