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