290 {
291 if (!event.
exists(ecal_sp_coll_name_, ecal_sp_pass_name_))
return;
293 ecal_sp_coll_name_, ecal_sp_pass_name_);
294
295
296 std::map<int, TruthParticle> by_id;
297 for (const auto& sp : sp_hits) {
298 auto pos = sp.getPosition();
299 auto mom = sp.getMomentum();
300 if (std::abs(240 - pos[2]) > 0.1 || mom[2] <= 0) continue;
301 double e = sp.getEnergy();
302 int apdg = std::abs(sp.getPdgID());
303 double vis = (apdg == 11 || apdg == 22)
304 ? e
305 : e - std::sqrt(std::max(0., e * e - mom[0] * mom[0] -
306 mom[1] * mom[1] -
307 mom[2] * mom[2]));
308 if (vis < truth_min_energy_) continue;
309 auto it = by_id.find(sp.getTrackID());
310 if (it != by_id.end() && it->second.visible_e_ >= vis) continue;
311 by_id[sp.getTrackID()] = {sp.getTrackID(), sp.getPdgID(), vis, e,
312 pos[0], pos[1], pos[2], mom[0],
313 mom[1], mom[2]};
314 }
315 std::vector<TruthParticle> truth;
316 double truth_e = 0;
317
318 for (auto& [id, t] : by_id) {
319 truth.push_back(t);
320 int cat = category(t.pdg_);
321 truth_e += (cat <= kP) ? t.total_e_ : t.visible_e_;
322 }
324 double pf_e = 0;
325 for (const auto& c : cands) pf_e += c.getEnergy();
326 if (truth_e > 0)
histograms_.
fill(
"sum_pf_e_over_truth", pf_e / truth_e);
327
328
329 std::vector<int> cand_truth(cands.size(), -1);
330 for (size_t i = 0; i < cands.size(); i++) {
331 if (!(cands[i].getPID() & 1)) continue;
332 for (size_t j = 0; j < truth.size(); j++) {
333 if (!truth[j].matched_ &&
334 truth[j].track_id_ == cands[i].getTruthTrackId()) {
335 truth[j].matched_ = true;
336 cand_truth[i] = j;
337 break;
338 }
339 }
340 }
341
342
343 std::vector<size_t> neutral;
344 for (size_t i = 0; i < cands.size(); i++) {
345 if (!(cands[i].getPID() & 1)) neutral.push_back(i);
346 }
347 std::sort(neutral.begin(), neutral.end(), [&cands](size_t a, size_t b) {
348 return cands[a].getEnergy() > cands[b].getEnergy();
349 });
350 for (size_t i : neutral) {
351 bool has_ecal = cands[i].getPID() & 2;
352 auto xyz =
353 has_ecal ? cands[i].getEcalClusterXYZ() : cands[i].getHcalClusterXYZ();
354
355
356 auto cand_dir = direction(cands[i]);
357 int best = -1;
358 double best_d = has_ecal ? neutral_match_dist_
359 : neutral_hcal_match_angle_ * M_PI / 180.;
360 for (size_t j = 0; j < truth.size(); j++) {
361 if (truth[j].matched_) continue;
362
363 double ratio = cands[i].getEnergy() / truth[j].visible_e_;
364 if (ratio < 1. / neutral_match_max_ratio_ ||
365 ratio > neutral_match_max_ratio_)
366 continue;
367 double d;
368 if (has_ecal) {
369 auto [tx, ty] = project(truth[j], xyz[2]);
370 d = std::hypot(tx - xyz[0], ty - xyz[1]);
371 } else {
372 const auto& t = truth[j];
373 double pt = std::hypot(t.px_, t.py_, t.pz_);
374 double c =
375 (t.px_ * cand_dir[0] + t.py_ * cand_dir[1] + t.pz_ * cand_dir[2]) /
376 pt;
377 d = std::acos(std::clamp(c, -1., 1.));
378 }
379 if (d < best_d) {
380 best_d = d;
381 best = j;
382 }
383 }
384 if (best >= 0) {
385 truth[best].matched_ = true;
386 cand_truth[i] = best;
387 }
388 }
389
390
391 std::vector<int> truth_cand(truth.size(), -1);
392 for (size_t i = 0; i < cands.size(); i++) {
393 int r_type = pfType(cands[i].getPdgId());
394 if (cand_truth[i] < 0) {
395
397 if (r_type == kPFH) {
399 massColumn(cands[i].getMass()));
400 }
401 if (!(cands[i].getPID() & 1)) {
403 }
404 continue;
405 }
406 truth_cand[cand_truth[i]] = i;
407 const auto& t = truth[cand_truth[i]];
408
409 double ref = (cands[i].getPID() & 1) ? t.total_e_ : t.visible_e_;
410 if (ref > 0) {
411 histograms_.
fill(
"response_vs_type", r_type, cands[i].getEnergy() / ref);
412 }
413 }
414
415
416 bool has_particles =
417 event.exists(sim_particles_coll_name_, sim_particles_pass_name_);
418 std::map<int, ldmx::SimParticle> particles;
419 if (has_particles) {
421 sim_particles_pass_name_);
422 }
423
424
425
426 std::vector<bool> seen(truth.size(), false);
427 for (size_t j = 0; j < truth.size(); j++) {
428 if (truth_cand[j] >= 0) continue;
429 for (const auto& cand : cands) {
430 int pid = cand.getPID();
431 if (pid & 4 && !(pid & 2)) {
432 auto d = direction(cand);
433 const auto& t = truth[j];
434 double pt = std::hypot(t.px_, t.py_, t.pz_);
435 double cs = (t.px_ * d[0] + t.py_ * d[1] + t.pz_ * d[2]) / pt;
436 if (cs > std::cos(neutral_hcal_match_angle_ * M_PI / 180.))
437 seen[j] = true;
438 } else {
439 auto xyz =
440 (pid & 2) ? cand.getEcalClusterXYZ() : cand.getEcalPositionXYZ();
441 auto [tx, ty] = project(truth[j], xyz[2]);
442 if (std::hypot(tx - xyz[0], ty - xyz[1]) < neutral_match_dist_)
443 seen[j] = true;
444 }
445 if (seen[j]) break;
446 }
447 }
448
449
450
451 for (size_t j = 0; j < truth.size(); j++) {
452 const auto& t = truth[j];
453 int c = truth_cand[j];
454 int r_type = c >= 0 ? pfType(cands[c].getPdgId())
455 : (seen[j] ? kPFOther : kPFNotFound);
456 bool reco = reconstructable(t.pdg_, t.px_, t.py_, t.pz_, t.visible_e_);
458 if (reco) {
460 r_type);
461 }
462 int sp = category(t.pdg_);
463 if (sp == kPi || sp == kK || sp == kP) {
464
465 int col = r_type == kPFH ? massColumn(cands[c].getMass())
466 : r_type == kPFNotFound ? 4
467 : 3;
469 if (reco) {
471 col);
472 }
473 }
474 if (sp == kNHad) {
475
476 int col = r_type == kPFNotFound ? 4 : 3;
477 if (r_type == kPFH0) col = 0;
478 if (r_type == kPFGamma) col = 1;
479 int region = kNoInteraction;
480 if (auto it = particles.find(t.track_id_); it != particles.end()) {
481 region = interactionRegion(it->second.getEndPoint());
482 }
483 int row = neutralRow(t.pdg_, region);
485 if (reco) {
487 }
488 }
489 }
490 for (const auto& t : truth) {
492 if (t.matched_)
histograms_.
fill(
"truth_type_matched", category(t.pdg_));
493 if (t.pdg_ == 22) {
495 if (t.matched_)
histograms_.
fill(
"truth_gamma_e_matched", t.visible_e_);
496 }
497 }
498
499
500 if (!has_particles) return;
501 auto pi0_parent = [&particles](int track_id) {
502 auto it = particles.find(track_id);
503 if (it == particles.end() || it->second.getParents().empty()) return -1;
504 int parent = it->second.getParents()[0];
505 auto pit = particles.find(parent);
506 if (pit == particles.end() || pit->second.getPdgID() != 111) return -1;
507 return parent;
508 };
509 for (size_t a = 0; a < cands.size(); a++) {
510 if (cands[a].getPdgId() != 22 || cand_truth[a] < 0) continue;
511 int pa = pi0_parent(truth[cand_truth[a]].track_id_);
512 if (pa < 0) continue;
513 for (size_t b = a + 1; b < cands.size(); b++) {
514 if (cands[b].getPdgId() != 22 || cand_truth[b] < 0) continue;
515 if (pi0_parent(truth[cand_truth[b]].track_id_) != pa) continue;
516 auto d1 = direction(cands[a]);
517 auto d2 = direction(cands[b]);
518 double cos12 = d1[0] * d2[0] + d1[1] * d2[1] + d1[2] * d2[2];
519 double m2 = 2 * cands[a].getEnergy() * cands[b].getEnergy() * (1 - cos12);
521 }
522 }
523
524 std::map<int, std::vector<int>> pi0_photons;
525 for (size_t j = 0; j < truth.size(); j++) {
526 if (truth[j].pdg_ != 22) continue;
527 int pa = pi0_parent(truth[j].track_id_);
528 if (pa >= 0) pi0_photons[pa].push_back(j);
529 }
530 for (const auto& [pa, photons] : pi0_photons) {
532 if (photons.size() < 2) continue;
533
534 int n_gamma = 0, n_h0 = 0, n_seen = 0;
535 for (int j : photons) {
536 int c = truth_cand[j];
537 n_seen += c >= 0 || seen[j];
538 if (c < 0) continue;
539 int r_type = pfType(cands[c].getPdgId());
540 n_gamma += r_type == kPFGamma;
541 n_h0 += r_type == kPFH0;
542 }
543 int col = n_seen > 0 ? 3 : 4;
544 if (n_gamma >= 2)
545 col = 2;
546 else if (n_h0 > 0)
547 col = 0;
548 else if (n_gamma == 1)
549 col = 1;
551 const auto& pi0 = particles.at(pa);
552 double pi0_ke = pi0.getEnergy() - pi0.getMass();
553 double pi0_theta = std::acos(
554 pi0.getMomentum()[2] /
555 std::max(1e-9, std::hypot(pi0.getMomentum()[0], pi0.getMomentum()[1],
556 pi0.getMomentum()[2])));
557 if (pi0_ke > reco_min_ke_pi0_ &&
558 pi0_theta < reco_max_theta_ * M_PI / 180.) {
560 }
561 }
562}
Class representing a simulated particle.
Represents a simulated tracker hit in the simulation.