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