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