Run the processor.
135 {
136 eventnr_++;
137
138 auto tg{geometry()};
139
140
141
142 std::vector<ldmx::Track> tracks;
143
144 auto start = std::chrono::high_resolution_clock::now();
145
146 nevents_++;
147
148 ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("LDMX Tracking Geometry Maker",
149 Acts::Logging::DEBUG));
150
151
152 Acts::PropagatorOptions<Acts::StepperPlainOptions,
153 Acts::NavigatorPlainOptions, ActionList>
154 propagator_options(geometryContext(), magneticFieldContext());
155
156 propagator_options.pathLimit = std::numeric_limits<double>::max();
157
158 propagator_options.loopProtection = false;
159
160
161
162 auto& m_interactor =
163 propagator_options.actorList.get<Acts::MaterialInteractor>();
164 m_interactor.multipleScattering = true;
165 m_interactor.energyLoss = true;
166 m_interactor.recordInteractions = false;
167
168
169 auto& s_logger =
170 propagator_options.actorList.get<Acts::detail::SteppingLogger>();
171 s_logger.sterile = true;
172
173 propagator_options.stepping.maxStepSize =
174 propagator_step_size_ * Acts::UnitConstants::mm;
175 propagator_options.maxSteps = propagator_max_steps_;
176
177
178
179
180
181
182
183
184
185 auto setup = std::chrono::high_resolution_clock::now();
186 profiling_map_["setup"] +=
187 std::chrono::duration<double, std::milli>(setup - start).count();
188
190 measurement_collection_, input_pass_name_);
191
192
193 std::shared_ptr<tracking::sim::TruthMatchingTool> truth_matching_tool =
194 nullptr;
195 std::map<int, ldmx::SimParticle> particle_map;
196
197 if (event.
exists(sim_particles_coll_name_, sim_particles_event_passname_)) {
198 ldmx_log(debug) << "Setting up track truth matching tool";
200 sim_particles_coll_name_, sim_particles_event_passname_);
201 truth_matching_tool = std::make_shared<tracking::sim::TruthMatchingTool>(
202 particle_map, measurements);
203 }
204
205
206
207 const auto geo_id_sl_map = makeGeoIdSourceLinkMap(tg, measurements);
208
209 auto hits = std::chrono::high_resolution_clock::now();
210 profiling_map_["hits"] +=
211 std::chrono::duration<double, std::milli>(hits - setup).count();
212
213
214
215
216 const auto& seed_tracks =
217 event.getCollection<
ldmx::Track>(seed_coll_name_, input_pass_name_);
218
219 ldmx_log(info) << "Number of " << seed_coll_name_
220 << " seed tracks = " << seed_tracks.size();
221
222 if (seed_tracks.empty()) {
223 std::vector<ldmx::Track> empty;
224 ldmx_log(warn) << "No seed tracks, returning...";
225 event.add(out_trk_collection_, empty);
226 return;
227 }
228
229
230 std::vector<Acts::BoundTrackParameters> start_parameters;
231
232 ldmx_log(debug) << "Transform the seed track to bound parameters";
233 int seed_track_index{0};
234 for (auto& seed : seed_tracks) {
235
236
237
238 Acts::Vector3 perigee_acts = tracking::sim::utils::ldmx2Acts(Acts::Vector3(
239 seed.getPerigeeX(), seed.getPerigeeY(), seed.getPerigeeZ()));
240 std::shared_ptr<Acts::PerigeeSurface> perigee_surface =
241 Acts::Surface::makeShared<Acts::PerigeeSurface>(perigee_acts);
242
243 Acts::BoundVector param_vec;
244 param_vec << seed.getD0(), seed.getZ0(), seed.getPhi(), seed.getTheta(),
245 seed.getQoP(), seed.getT();
246
247 Acts::BoundMatrix cov_mat =
248 tracking::sim::utils::unpackCov(seed.getPerigeeCov());
249
250 ldmx_log(debug) << " For seed index_ = " << seed_track_index
251 << ": Perigee X / Y / Z = " << seed.getPerigeeX() << " / "
252 << seed.getPerigeeY() << " / " << seed.getPerigeeZ()
253 << ", D0 = " << param_vec[0] << ", Z0 = " << param_vec[1]
254 << ", Phi = " << param_vec[2]
255 << ", Theta = " << param_vec[3]
256 << ", QoP = " << param_vec[4]
257 << ", Time = " << param_vec[5];
258
259 ldmx_log(debug) << " Cov matrix diagonal (" << cov_mat(0, 0) << ", "
260 << cov_mat(1, 1) << ", " << cov_mat(2, 2) << ")";
261
262
263 auto part_hypo{Acts::ParticleHypothesis::electron()};
264 start_parameters.push_back(Acts::BoundTrackParameters(
265 perigee_surface, param_vec, cov_mat, part_hypo));
266
267
269
270 seed_track_index++;
271 }
272
273 auto seeds = std::chrono::high_resolution_clock::now();
274 profiling_map_["seeds"] +=
275 std::chrono::duration<double, std::milli>(seeds - hits).count();
276
277 Acts::GainMatrixUpdater kf_updater;
278
279
280
281
282 Acts::MeasurementSelector::Config measurement_selector_cfg = {
283
284 {Acts::GeometryIdentifier(), {{}, {outlier_pval_}, {1u}}},
285 };
286
287 Acts::MeasurementSelector meas_sel{measurement_selector_cfg};
288
290
291
292 struct SourceLinkAccIt {
293 using BaseIt = decltype(geo_id_sl_map.begin());
294 BaseIt it_;
295
296#pragma GCC diagnostic push
297#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
298
299 using difference_type = typename BaseIt::difference_type;
300 using iterator_category = typename BaseIt::iterator_category;
301 using value_type = Acts::SourceLink;
302 using pointer = typename BaseIt::pointer;
303 using reference = value_type&;
304#pragma GCC diagnostic pop
305
306 SourceLinkAccIt& operator++() {
307 ++it_;
308 return *this;
309 }
310 bool operator==(const SourceLinkAccIt& other) const {
311 return it_ == other.it_;
312 }
313 bool operator!=(const SourceLinkAccIt& other) const {
314 return !(*this == other);
315 }
316 value_type operator*() const { return value_type{it_->second}; }
317 };
318
319 auto source_link_accessor = [&](const Acts::Surface& surface)
320 -> std::pair<SourceLinkAccIt, SourceLinkAccIt> {
321 auto [begin, end] = geo_id_sl_map.equal_range(surface.geometryId());
322 return {SourceLinkAccIt{begin}, SourceLinkAccIt{end}};
323 };
324
325
326 Acts::TrackStateCreator<SourceLinkAccIt, TrackContainer> track_state_creator;
327 track_state_creator.sourceLinkAccessor
328 .connect<&decltype(source_link_accessor)::operator(),
329 decltype(source_link_accessor)>(&source_link_accessor);
330 if (use1_dmeasurements_) {
331 track_state_creator.calibrator
333 Acts::VectorMultiTrajectory>>(&calibrator);
334 } else {
335 track_state_creator.calibrator
337 Acts::VectorMultiTrajectory>>(&calibrator);
338 }
339 track_state_creator.measurementSelector
340 .connect<&Acts::MeasurementSelector::select<Acts::VectorMultiTrajectory>>(
341 &meas_sel);
342
343 Acts::CombinatorialKalmanFilterExtensions<TrackContainer> ckf_extensions;
344 ckf_extensions.updater.connect<
345 &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(
346 &kf_updater);
347 ckf_extensions.createTrackStates.connect<&Acts::TrackStateCreator<
348 SourceLinkAccIt, TrackContainer>::createTrackStates>(
349 &track_state_creator);
350
351 ldmx_log(debug) << "Setting up surfaces...";
352
353 std::shared_ptr<const Acts::PerigeeSurface> origin_surface =
354 Acts::Surface::makeShared<Acts::PerigeeSurface>(
355 Acts::Vector3(0., 0., 0.));
356
357 ldmx_log(debug) << "About to run CKF...";
358
359
360 auto ckf_setup = std::chrono::high_resolution_clock::now();
361 profiling_map_["ckf_setup"] +=
362 std::chrono::duration<double, std::milli>(ckf_setup - seeds).count();
363
364 Acts::VectorTrackContainer vtc;
365 Acts::VectorMultiTrajectory mtj;
366 Acts::TrackContainer tc{vtc, mtj};
367
368
369
370 ldmx_log(debug) << "Loop on the track candidates";
371 for (size_t track_id = 0u; track_id < start_parameters.size(); ++track_id) {
372 ldmx_log(debug) << "---------------------------";
373 ldmx_log(debug) << "Candidate Track ID = " << track_id;
374
375 const Acts::CombinatorialKalmanFilterOptions<TrackContainer> ckf_options(
376 TrackingGeometryUser::geometryContext(),
377 TrackingGeometryUser::magneticFieldContext(),
378 TrackingGeometryUser::calibrationContext(), ckf_extensions,
379 static_cast<Acts::PropagatorPlainOptions>(propagator_options),
380 true , false );
381
382 ldmx_log(debug) << " Checking options: multiple scattering = "
383 << ckf_options.multipleScattering
384 << " energy loss = " << ckf_options.energyLoss;
385
386
387 auto results =
388 ckf_->findTracks(start_parameters.at(track_id), ckf_options, tc);
389
390 auto start_params = start_parameters.at(track_id).parameters().transpose();
391
392
393 if (!results.ok()) {
394 if (!tagger_tracking_) {
395
396 n_fieldmap_ckf_failed_recoil_++;
397 ldmx_log(debug)
398 << " Field-map CKF failed, trying zero-B CKF fallback";
399 results = ckf_zero_b_->findTracks(start_parameters.at(track_id),
400 ckf_options, tc);
401 if (results.ok()) {
402 n_zerob_ckf_recovered_recoil_++;
403 ldmx_log(debug) << " Yay! Zero-B CKF succeeded as fallback!";
404 } else {
405 ldmx_log(debug) << " Zero-B CKF also failed!";
406 }
407 } else {
408
409 n_fieldmap_ckf_failed_tagger_++;
410 ldmx_log(debug)
411 << " Field-map CKF failed, trying const-B (1.5T) CKF fallback";
412 results = ckf_const_b_->findTracks(start_parameters.at(track_id),
413 ckf_options, tc);
414 if (results.ok()) {
415 n_constb_ckf_recovered_tagger_++;
416 ldmx_log(debug) << " Yay! Const-B CKF succeeded as fallback!";
417 } else {
418 ldmx_log(debug) << " Const-B CKF also failed!";
419 }
420 }
421 }
422
423 ldmx_log(debug)
424 << " Checking CKF success for track candidate with params: "
425 << " D0 = " << start_params[0] << " Z0 = " << start_params[1]
426 << ", Phi = " << start_params[2] << " Theta = " << start_params[3]
427 << ", QoP = " << start_params[4] << " Time = " << start_params[5];
428 if (not results.ok()) {
429 ldmx_log(debug) << " CKF failed!";
430 continue;
431 } else {
432 ldmx_log(debug) << " CKF succeded!";
433 }
434
435 auto& tracks_from_seed = results.value();
436 if (tracks_from_seed.size() != 1) {
437 ldmx_log(info) << " tracksFromSeed.size = " << tracks_from_seed.size();
438 }
439
440 for (auto& track : tracks_from_seed) {
441
442 auto smooth_result = Acts::smoothTrack(geometryContext(), track);
443 if (!smooth_result.ok()) {
444 ldmx_log(warn) << "smoothTrack failed: "
445 << smooth_result.error().message();
446 }
447
449
450
451 auto opt_target = trk_extrap_->extrapolate(track, target_surface_);
452
453 if (!opt_target) {
454 if (tagger_tracking_) {
455 n_fieldmap_target_extrap_failed_tagger_++;
456 ldmx_log(debug) << " Field-map target extrapolation failed, "
457 "trying const-B (1.5T) fallback";
458 opt_target = trk_extrap_const_b_->extrapolate(track, target_surface_);
459 if (opt_target)
460 n_constb_target_extrap_recovered_tagger_++;
461 else
462 ldmx_log(debug) << " Both field-map and Const-B target "
463 "extrapolation failed!";
464 } else {
465 n_fieldmap_target_extrap_failed_recoil_++;
466 ldmx_log(debug) << " Field-map target extrapolation failed, "
467 "trying zero-B fallback";
468 opt_target = trk_extrap_zero_b_->extrapolate(track, target_surface_);
469 if (opt_target)
470 n_zerob_target_extrap_recovered_recoil_++;
471 else
472 ldmx_log(debug)
473 << " Both field-map and Zero-B target extrapolation failed!";
474 }
475 }
476
477 if (!opt_target) {
478 ldmx_log(debug) << " Could not extrapolate to target! nhits = "
479 << track.nMeasurements() << " Printing track states:";
480 for (const auto ts : track.trackStatesReversed()) {
481 if (ts.hasSmoothed())
482 ldmx_log(debug) << " Parameters: " << ts.smoothed().transpose();
483 else
484 ldmx_log(debug) << " Track state not smoothed!";
485 }
486 ldmx_log(debug) << " ...skipping this track candidate...";
487 continue;
488 }
489
490 ldmx_log(debug) << " Successfully obtained TrackState at target";
491
492
493 auto ts_at_target = tracking::sim::utils::makeTrackState(
494 geometryContext(), *opt_target, ldmx::AtTarget);
495 trk.addTrackState(ts_at_target);
496
497 ldmx_log(debug) << " Position at target (LDMX): ("
498 << ts_at_target.pos_[0] << ", " << ts_at_target.pos_[1]
499 << ", " << ts_at_target.pos_[2] << ") mm"
500 << " Momentum: (" << ts_at_target.mom_[0] << ", "
501 << ts_at_target.mom_[1] << ", " << ts_at_target.mom_[2]
502 << ") GeV";
503
504
505 track.setReferenceSurface(target_surface_);
506 track.parameters() = opt_target->parameters();
507
508
509 trk.setPerigeeParameters(tracking::sim::utils::convertActsToLdmxPars(
510 opt_target->parameters()));
511 if (opt_target->covariance()) {
512 std::vector<double> cov_vec;
513 tracking::sim::utils::flatCov(*(opt_target->covariance()), cov_vec);
514 trk.setPerigeeCov(cov_vec);
515 }
516
517 Acts::Vector3 target_loc_ldmx = tracking::sim::utils::acts2Ldmx(
518 target_surface_->localToGlobalTransform(geometryContext())
519 .translation());
520 trk.setPerigeeLocation(target_loc_ldmx[0], target_loc_ldmx[1],
521 target_loc_ldmx[2]);
522
523 trk.setChi2(track.chi2());
524 trk.setNhits(track.nMeasurements());
525 trk.setNdf(track.nMeasurements() - 5);
526 trk.setNsharedHits(track.nSharedHits());
527 trk.setCharge(opt_target->parameters()[Acts::eBoundQOverP] > 0 ? 1 : -1);
528
529
530 if ((trk.getNhits() <= min_hits_) ||
531 (std::abs(1. / trk.getQoP()) <= 0.05)) {
532 ldmx_log(debug)
533 << " > Track candidate did NOT meet the requirements: Nhits = "
534 << trk.getNhits() << " and p = " << std::abs(1. / trk.getQoP())
535 << " GeV";
536 continue;
537 }
538
539
540 ldmx_log(debug) << " Add measurements to the final track from "
541 << track.nTrackStates() << " TrackStates with "
542 << track.nMeasurements() << " measurements";
543
544 int trk_state_index{0};
545 for (const auto ts : track.trackStatesReversed()) {
546
547 ldmx_log(debug) << " Checking Track State index_ = "
548 << trk_state_index << " at location "
549 << ts.referenceSurface()
550 .localToGlobalTransform(geometryContext())
551 .translation()
552 .transpose();
553
554 if (ts.hasSmoothed()) {
555 ldmx_log(debug) << " Smoothed track parameters: "
556 << ts.smoothed().transpose();
557
558
559 }
560
561
562 auto type_flags = ts.typeFlags();
563
564 if (type_flags.isMeasurement() && ts.hasUncalibratedSourceLink()) {
565 Acts::SourceLink usl = ts.getUncalibratedSourceLink();
568
570 ldmx_log(debug) << " Adding measurement to ldmx::track with "
571 "source link index_ = "
573 ldmx_log(trace) << " Measurement:\n" << ldmx_meas;
574 trk.addMeasurementIndex(sl.
index());
575
576
577
578
579
580
581
582 if (ts.hasSmoothed()) {
583 trk.addSmoothedLoc0(
584 static_cast<float>(ts.smoothed()[Acts::eBoundLoc0]),
585 static_cast<float>(ts.smoothedCovariance()(Acts::eBoundLoc0,
586 Acts::eBoundLoc0)));
587 }
588
589
590 if (ts.hasSmoothed()) {
591 const auto& meas_surface = ts.referenceSurface();
592 const auto& smoothed_params = ts.smoothed();
593
594
595
596
597 float p_inv = smoothed_params[Acts::eBoundQOverP];
598 float p = 1.0f / std::abs(p_inv);
599 float theta = smoothed_params[Acts::eBoundTheta];
600 float phi = smoothed_params[Acts::eBoundPhi];
601
602 Acts::Vector3 global_momentum(p * std::sin(theta) * std::cos(phi),
603 p * std::sin(theta) * std::sin(phi),
604 p * std::cos(theta));
605
606
607 auto local_frame_transform =
608 meas_surface.localToGlobalTransform(geometryContext());
609 Acts::Vector3 local_momentum =
610 local_frame_transform.rotation().transpose() * global_momentum;
611
612
613 float phi_u = (local_momentum.z() != 0)
614 ? local_momentum.x() / local_momentum.z()
615 : 0.;
616 float phi_v = (local_momentum.z() != 0)
617 ? local_momentum.y() / local_momentum.z()
618 : 0.;
619
620
621
622
623
624 float sensor_thickness = 0.0f;
625 if (const auto* placement = meas_surface.surfacePlacement()) {
626 sensor_thickness = static_cast<float>(
628 ->thickness());
629 } else {
630 ldmx_log(warn) << "No detector element for measurement surface"
631 << " — skipping dE/dx for this hit";
632 continue;
633 }
634 float tan_angle_sq = phi_u * phi_u + phi_v * phi_v;
635 float cos_angle = 1.0f / std::sqrt(1.0f + tan_angle_sq);
636 float path_length = sensor_thickness / cos_angle;
637
638 ldmx_log(debug) << " Local angles: phi_u = " << phi_u
639 << ", phi_v = " << phi_v
640 << "; Path length = " << path_length << " mm";
641
642
643 float edep = ldmx_meas.getEdep();
644 float dedx = edep / path_length;
645 trk.addDedxMeasurement(dedx);
646
647 ldmx_log(debug) << " Edep = " << edep
648 << " MeV, dE/dx = " << dedx << " MeV/mm";
649 }
650 } else {
651 ldmx_log(debug) << " This TrackState is not a measurement";
652 }
653 trk_state_index++;
654 }
655
656 ldmx_log(debug) << " Starting extrapolations";
657
658
659 const double ecal_scoring_plane = 240.5;
660 Acts::Vector3 pos(ecal_scoring_plane, 0., 0.);
661 Acts::Translation3 surf_translation(pos);
662 Acts::Transform3 surf_transform(surf_translation * surf_rotation_);
663 const std::shared_ptr<Acts::PlaneSurface> ecal_surface =
664 Acts::Surface::makeShared<Acts::PlaneSurface>(surf_transform);
665
666
667 const std::shared_ptr<Acts::Surface> beam_origin_surface =
668 tracking::sim::utils::unboundSurface(-700);
669
670 if (tagger_tracking_) {
671 ldmx_log(debug) << " Beam Origin Extrapolation";
672 auto opt_beam_origin =
673 trk_extrap_->extrapolate(track, beam_origin_surface);
674 if (opt_beam_origin) {
675 trk.addTrackState(tracking::sim::utils::makeTrackState(
676 geometryContext(), *opt_beam_origin, ldmx::AtBeamOrigin));
677 ldmx_log(debug)
678 << " Successfully obtained TrackState at beam origin";
679 }
680 }
681
682
683 if (!tagger_tracking_) {
684 ldmx_log(debug) << " Ecal Extrapolation";
685 auto opt_ecal = trk_extrap_->extrapolate(track, ecal_surface);
686
687 if (!opt_ecal) {
688 n_fieldmap_ecal_extrap_failed_recoil_++;
689 ldmx_log(debug) << " Field-map ECAL extrapolation failed, trying "
690 "zero-B fallback";
691 opt_ecal = trk_extrap_zero_b_->extrapolate(track, ecal_surface);
692 if (opt_ecal)
693 n_zerob_ecal_extrap_recovered_recoil_++;
694 else
695 ldmx_log(debug)
696 << " Both field-map and Zero-B ECAL extrapolation failed!";
697 }
698
699 if (opt_ecal) {
700 auto ts_at_ecal = tracking::sim::utils::makeTrackState(
701 geometryContext(), *opt_ecal, ldmx::AtECAL);
702 trk.addTrackState(ts_at_ecal);
703 ldmx_log(debug) << " Successfully obtained TrackState at ECAL";
704 ldmx_log(debug) << " Position at ECAL (LDMX): ("
705 << ts_at_ecal.pos_[0] << ", " << ts_at_ecal.pos_[1]
706 << ", " << ts_at_ecal.pos_[2] << ") mm";
707 }
708 }
709
710
711 if (truth_matching_tool) {
712 auto truth_info = truth_matching_tool->truthMatch(trk);
713 trk.setTrackID(truth_info.track_id_);
714 trk.setPdgID(truth_info.pdg_id_);
715 trk.setTruthProb(truth_info.truth_prob_);
716 }
717
718
719 ldmx_log(debug)
720 << " > Adding the track candidate to the track collection";
721 tracks.push_back(trk);
722 ntracks_++;
723 }
724 }
725
726 ldmx_log(info) << "Number of CKF tracks " << tracks.size();
727
728 auto ckf_run = std::chrono::high_resolution_clock::now();
729 profiling_map_["ckf_run"] +=
730 std::chrono::duration<double, std::milli>(ckf_run - ckf_setup).count();
731
732
733 auto shared_hits = computeSharedHits(
734 tracks, measurements, tg, tracking::sim::utils::sourceLinkHash,
735 tracking::sim::utils::sourceLinkEquality);
736 for (std::size_t i_track = 0; i_track < shared_hits.size(); ++i_track) {
737 tracks[i_track].setNsharedHits(shared_hits[i_track].size());
738 for (auto idx : shared_hits[i_track]) {
739 tracks[i_track].addSharedIndex(idx);
740 }
741 }
742
743 auto result_loop = std::chrono::high_resolution_clock::now();
744 profiling_map_["result_loop"] +=
745 std::chrono::duration<double, std::milli>(result_loop - ckf_run).count();
746
747
748 event.add(out_trk_collection_, tracks);
749
750 auto end = std::chrono::high_resolution_clock::now();
751
752
753 auto diff = end - start;
754 processing_time_ += std::chrono::duration<double, std::milli>(diff).count();
755}
constexpr Index index() const
Access the index_.
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.
Class representing a simulated particle.
Implementation of a track object.
void calibrate1d(const Acts::GeometryContext &, const Acts::CalibrationContext &, const Acts::SourceLink &genericSourceLink, typename traj_t::TrackStateProxy trackState) const
Find the measurement corresponding to the source link.
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.