Process event to find ECAL tracks.
521 {
522 auto start = std::chrono::high_resolution_clock::now();
523 nevents_++;
524
525 std::vector<ldmx::Track> tracks;
526
527
528 if (!event.
exists(rec_coll_name_, rec_pass_name_)) {
529 ldmx_log(debug) << "No ECAL RecHits collection found";
530 event.add(out_track_collection_, tracks);
531 return;
532 }
533
534 const std::vector<ldmx::EcalHit> ecal_hits =
535 event.getCollection<
ldmx::EcalHit>(rec_coll_name_, rec_pass_name_);
536
537 ldmx_log(debug) << "Processing " << ecal_hits.size() << " ECAL hits";
538
539
540 std::vector<double> measurement_energies;
542 ldmx_log(debug) << "Created " << measurements.size() << " measurements";
543
544 if (measurements.empty()) {
545 event.add(out_track_collection_, tracks);
546 return;
547 }
548
549
551 ldmx_log(info) << "Source link map: " << geo_id_sl_map.size()
552 << " entries from " << measurements.size() << " measurements";
553
554
555 auto seed_tracks =
findSeeds(measurements);
556 ldmx_log(info) << "Found " << seed_tracks.size() << " seed tracks";
557 nseeds_ += seed_tracks.size();
558
559 if (seed_tracks.empty()) {
560 event.add(out_track_collection_, tracks);
561 return;
562 }
563
564
567 .get();
570 .get();
573 .get();
574
575
576 Acts::PropagatorPlainOptions propagator_options(gctx, mctx);
577 propagator_options.pathLimit = std::numeric_limits<double>::max();
578 propagator_options.maxSteps = 1000;
579 propagator_options.stepping.maxStepSize = 100.0 * Acts::UnitConstants::mm;
580
581
582 Acts::GainMatrixUpdater kf_updater;
583 Acts::MeasurementSelector::Config meas_sel_cfg = {
584 {Acts::GeometryIdentifier(), {{}, {max_chi2_}, {1u}}}};
585 Acts::MeasurementSelector meas_sel{meas_sel_cfg};
586
588
589
590 struct SourceLinkAccIt {
591 using BaseIt = decltype(geo_id_sl_map.begin());
592 BaseIt it_;
593
594#pragma GCC diagnostic push
595#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
596
597 using difference_type = typename BaseIt::difference_type;
598 using iterator_category = std::input_iterator_tag;
599 using value_type = Acts::SourceLink;
600 using pointer = value_type*;
601 using reference = value_type&;
602#pragma GCC diagnostic pop
603
604 SourceLinkAccIt& operator++() {
605 ++it_;
606 return *this;
607 }
608 bool operator==(const SourceLinkAccIt& other) const {
609 return it_ == other.it_;
610 }
611 bool operator!=(const SourceLinkAccIt& other) const {
612 return !(*this == other);
613 }
614 value_type operator*() const { return value_type{it_->second}; }
615 };
616
617 auto source_link_accessor = [&](const Acts::Surface& surface)
618 -> std::pair<SourceLinkAccIt, SourceLinkAccIt> {
619 auto [begin, end] = geo_id_sl_map.equal_range(surface.geometryId());
620 return {SourceLinkAccIt{begin}, SourceLinkAccIt{end}};
621 };
622
623
624 Acts::TrackStateCreator<SourceLinkAccIt, TrackContainer> track_state_creator;
625 track_state_creator.sourceLinkAccessor
626 .connect<&decltype(source_link_accessor)::operator(),
627 decltype(source_link_accessor)>(&source_link_accessor);
628 track_state_creator.calibrator
630 Acts::VectorMultiTrajectory>>(&calibrator);
631 track_state_creator.measurementSelector
632 .connect<&Acts::MeasurementSelector::select<Acts::VectorMultiTrajectory>>(
633 &meas_sel);
634
635 Acts::CombinatorialKalmanFilterExtensions<TrackContainer> ckf_extensions;
636 ckf_extensions.updater.connect<
637 &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(
638 &kf_updater);
639 ckf_extensions.createTrackStates.connect<&Acts::TrackStateCreator<
640 SourceLinkAccIt, TrackContainer>::createTrackStates>(
641 &track_state_creator);
642
643
644 Acts::VectorTrackContainer vtc;
645 Acts::VectorMultiTrajectory mtj;
646 Acts::TrackContainer tc{vtc, mtj};
647
648
649 for (size_t seed_idx = 0; seed_idx < seed_tracks.size(); ++seed_idx) {
650 const auto& seed = seed_tracks[seed_idx];
651
652
653
654
655
656 Acts::BoundVector param_vec;
657 param_vec << seed.getD0(), seed.getZ0(), seed.getPhi(), seed.getTheta(),
658 seed.getQoP(), seed.getT();
659
660 ldmx_log(debug) << "Seed " << seed_idx << ": loc0=" << param_vec[0]
661 << " loc1=" << param_vec[1] << " phi=" << param_vec[2]
662 << " theta=" << param_vec[3] << " qop=" << param_vec[4];
663
664 Acts::BoundMatrix cov_mat =
665 tracking::sim::utils::unpackCov(seed.getPerigeeCov());
666
667 auto part_hypo{Acts::ParticleHypothesis::electron()};
668 auto& layer0_surface = layer_surfaces_.begin()->second;
669 Acts::BoundTrackParameters start_params(layer0_surface, param_vec, cov_mat,
670 part_hypo);
671
672
673 const Acts::CombinatorialKalmanFilterOptions<TrackContainer> ckf_options(
674 gctx, mctx, cctx, ckf_extensions, propagator_options);
675
676
677 auto results = ckf_->findTracks(start_params, ckf_options, tc);
678
679 if (!results.ok()) {
680 ldmx_log(debug) << "CKF failed for seed " << seed_idx << ": "
681 << results.error().message();
682 continue;
683 }
684
685 auto& tracks_from_seed = results.value();
686 ldmx_log(info) << "CKF returned " << tracks_from_seed.size()
687 << " tracks from seed " << seed_idx;
688 for (auto& track : tracks_from_seed) {
689
690 int n_meas = 0, n_holes = 0, n_outliers = 0, n_total = 0;
691 for (const auto& ts : track.trackStatesReversed()) {
692 ++n_total;
693 if (ts.typeFlags().isMeasurement()) ++n_meas;
694 if (ts.typeFlags().isHole()) ++n_holes;
695 if (ts.typeFlags().isOutlier()) ++n_outliers;
696 }
697 ldmx_log(info) << "Track states: total=" << n_total << " meas=" << n_meas
698 << " holes=" << n_holes << " outliers=" << n_outliers;
699
700
701 auto smooth_result = Acts::smoothTrack(gctx, track);
702 if (!smooth_result.ok()) {
703 ldmx_log(warn) << "smoothTrack failed: "
704 << smooth_result.error().message();
705 continue;
706 }
707
708
710
711
712
713
714 Acts::BoundVector smoothed_params;
715 std::shared_ptr<const Acts::Surface> smoothed_surface;
716 bool found_smoothed = false;
717
718 for (const auto& ts : track.trackStatesReversed()) {
719 if (ts.hasSmoothed()) {
720 smoothed_params = ts.smoothed();
721 smoothed_surface = ts.referenceSurface().getSharedPtr();
722 found_smoothed = true;
723 break;
724 }
725 }
726
727 if (!found_smoothed) {
728 ldmx_log(warn) << "No smoothed track state found after smoothing";
729 continue;
730 }
731
732 ldmx_log(debug) << "Smoothed params: loc0=" << smoothed_params[0]
733 << " loc1=" << smoothed_params[1]
734 << " phi=" << smoothed_params[2]
735 << " theta=" << smoothed_params[3]
736 << " qop=" << smoothed_params[4];
737
738
739 Acts::FreeVector free_params = Acts::transformBoundToFreeParameters(
740 *smoothed_surface, gctx, smoothed_params);
741
742
743 Acts::Vector3 pos_acts(free_params[Acts::eFreePos0],
744 free_params[Acts::eFreePos1],
745 free_params[Acts::eFreePos2]);
746 Acts::Vector3 mom_acts(free_params[Acts::eFreeDir0],
747 free_params[Acts::eFreeDir1],
748 free_params[Acts::eFreeDir2]);
749
750 Acts::Vector3 pos_ldmx = tracking::sim::utils::acts2Ldmx(pos_acts);
751 Acts::Vector3 mom_ldmx = tracking::sim::utils::acts2Ldmx(mom_acts);
752
753 double x = pos_ldmx[0];
754 double y = pos_ldmx[1];
755 double z = pos_ldmx[2];
756 double px = mom_ldmx[0];
757 double py = mom_ldmx[1];
758 double pz = mom_ldmx[2];
759
760 ldmx_log(debug) << "LDMX momentum: px=" << px << " py=" << py
761 << " pz=" << pz;
762
763
764 double pt = std::sqrt(px * px + py * py);
765 double theta = std::atan2(pt, pz);
766 double phi = std::atan2(py, px);
767 if (phi < 0)
768 phi += 2.0 * M_PI;
769 double qop = free_params[Acts::eFreeQOverP];
770
771
772
773
774 double d0 = -(x * std::sin(phi) - y * std::cos(phi));
775 double z0 = z;
776 double time = free_params[Acts::eFreeTime];
777
778 Acts::BoundVector perigee_params;
779 perigee_params << d0, z0, phi, theta, qop, time;
780
781 ldmx_log(debug) << "Perigee params: d0=" << d0 << " z0=" << z0
782 << " phi=" << phi << " theta=" << theta;
783
784
785 trk.setPerigeeParameters(
786 tracking::sim::utils::convertActsToLdmxPars(perigee_params));
787
788
789 std::vector<double> cov_vec;
790 tracking::sim::utils::flatCov(track.covariance(), cov_vec);
791 trk.setPerigeeCov(cov_vec);
792
793 Acts::Vector3 ref_loc_ldmx = tracking::sim::utils::acts2Ldmx(
794 layer_surfaces_.begin()->second->center(gctx));
795 trk.setPerigeeLocation(ref_loc_ldmx[0], ref_loc_ldmx[1], ref_loc_ldmx[2]);
796
797 trk.setChi2(track.chi2());
798 trk.setNhits(track.nMeasurements());
799 trk.setNdf(track.nMeasurements() - 5);
800 trk.setNsharedHits(0);
801 trk.setCharge(qop > 0 ? 1 : -1);
802
803
804 for (const auto ts : track.trackStatesReversed()) {
805 if (ts.typeFlags().isMeasurement() && ts.hasUncalibratedSourceLink()) {
806 Acts::SourceLink usl = ts.getUncalibratedSourceLink();
809 trk.addMeasurementIndex(sl.
index());
810 }
811 }
812
813
814 double track_energy = 0.0;
815
816 if (use_roc_energy_ && !roc_range_values_.empty()) {
817
818
819
820
821 double trk_p_mag = 1.0 / std::abs(smoothed_params[Acts::eBoundQOverP]);
822 double trk_theta_deg = theta * 180.0 / M_PI;
823
824
825 std::vector<float> ele_radii(roc_range_values_[0].begin() + 4,
826 roc_range_values_[0].end());
827 for (const auto& row : roc_range_values_) {
828 float theta_min = row[0], theta_max = row[1];
829 float p_min = row[2], p_max = row[3];
830 bool inrange = true;
831 if (theta_min != -1.0f)
832 inrange = inrange && (trk_theta_deg >= theta_min);
833 if (theta_max != -1.0f)
834 inrange = inrange && (trk_theta_deg < theta_max);
835 if (p_min != -1.0f) inrange = inrange && (trk_p_mag >= p_min);
836 if (p_max != -1.0f) inrange = inrange && (trk_p_mag < p_max);
837 if (inrange) {
838 ele_radii.assign(row.begin() + 4, row.end());
839 }
840 }
841
842
843
844 for (const auto& hit : ecal_hits) {
845 if (hit.isNoise()) continue;
847 int layer = ecal_id.layer();
848 if (layer < 0 || layer >= static_cast<int>(ele_radii.size()))
849 continue;
850
851 auto [hx, hy, hz] = geometry_->
getPosition(ecal_id);
852
853
854 double dz = hz - z;
855 double proj_x = x + (px / pz) * dz;
856 double proj_y = y + (py / pz) * dz;
857
858
859 double dx = hx - proj_x;
860 double dy = hy - proj_y;
861 double dist = std::sqrt(dx * dx + dy * dy);
862
863 if (dist < ele_radii[layer]) {
864 track_energy += hit.getEnergy();
865 }
866 }
867 } else {
868
869 for (const auto ts : track.trackStatesReversed()) {
870 if (ts.typeFlags().isMeasurement() &&
871 ts.hasUncalibratedSourceLink()) {
872 Acts::SourceLink usl = ts.getUncalibratedSourceLink();
875 track_energy += measurement_energies[sl.
index()];
876 }
877 }
878 }
879
880
881 double track_p = track_energy;
882 qop = (track_p > 0) ? -1.0 / track_p : 0.0;
883 perigee_params[Acts::eBoundQOverP] = qop;
884 trk.setPerigeeParameters(
885 tracking::sim::utils::convertActsToLdmxPars(perigee_params));
886
887 ldmx_log(debug) << "Track energy from RecHits: " << track_energy
888 << " MeV, q/p=" << qop;
889
890 tracks.push_back(trk);
891 ntracks_++;
892 }
893 }
894
895 ldmx_log(info) << "Found " << tracks.size() << " fitted tracks";
896
897
898 event.add(out_track_collection_, tracks);
899
900 auto end = std::chrono::high_resolution_clock::now();
901 auto diff = end - start;
902 processing_time_ += std::chrono::duration<double, std::milli>(diff).count();
903}
constexpr Index index() const
Access the index_.
std::vector< ldmx::Measurement > createMeasurements(const std::vector< ldmx::EcalHit > &hits, std::vector< double > &energies)
Create ACTS measurement objects from ECAL hits.
std::unordered_multimap< Acts::GeometryIdentifier, acts_examples::IndexSourceLink > makeGeoIdSourceLinkMap(const std::vector< ldmx::Measurement > &measurements)
Create geometry ID to source link map for CKF.
std::vector< ldmx::Track > findSeeds(const std::vector< ldmx::Measurement > &measurements)
Find seed tracks via straight-line fitting.
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.
Stores reconstructed hit information from the ECAL.
static const std::string NAME
Conditions object name.
static const std::string NAME
Conditions object name.
void calibrate(const Acts::GeometryContext &, const Acts::CalibrationContext &, const Acts::SourceLink &genericSourceLink, typename traj_t::TrackStateProxy trackState) const
Find the measurement corresponding to the source link.