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