Run the processor.
152 {
153 eventnr_++;
154
155 auto tg{geometry()};
156
157
158
159 std::vector<ldmx::Track> tracks;
160
161 auto start = std::chrono::high_resolution_clock::now();
162
163 nevents_++;
164
165 auto logging_level = Acts::Logging::DEBUG;
166 ACTS_LOCAL_LOGGER(
167 Acts::getDefaultLogger("LDMX Tracking Geometry Maker", logging_level));
168
169
170 Acts::PropagatorOptions<Acts::StepperPlainOptions,
171 Acts::NavigatorPlainOptions, ActionList, AbortList>
172 propagator_options(geometryContext(), magneticFieldContext());
173
174 propagator_options.pathLimit = std::numeric_limits<double>::max();
175
176 propagator_options.loopProtection = false;
177
178
179
180 auto& m_interactor =
181 propagator_options.actionList.get<Acts::MaterialInteractor>();
182 m_interactor.multipleScattering = true;
183 m_interactor.energyLoss = true;
184 m_interactor.recordInteractions = false;
185
186
187 auto& s_logger =
188 propagator_options.actionList.get<Acts::detail::SteppingLogger>();
189 s_logger.sterile = true;
190
191 propagator_options.stepping.maxStepSize =
192 propagator_step_size_ * Acts::UnitConstants::mm;
193 propagator_options.maxSteps = propagator_max_steps_;
194
195
196
197
198
199
200
201
202
203 auto setup = std::chrono::high_resolution_clock::now();
204 profiling_map_["setup"] +=
205 std::chrono::duration<double, std::milli>(setup - start).count();
206
207 const std::vector<ldmx::Measurement> measurements =
209 input_pass_name_);
210
211
212 std::shared_ptr<tracking::sim::TruthMatchingTool> truth_matching_tool =
213 nullptr;
214 std::map<int, ldmx::SimParticle> particle_map;
215
216 if (event.
exists(
"SimParticles", sim_particles_event_passname_)) {
217 ldmx_log(debug) << "Setting up track truth matching tool";
219 "SimParticles", sim_particles_event_passname_);
220 truth_matching_tool = std::make_shared<tracking::sim::TruthMatchingTool>(
221 particle_map, measurements);
222 }
223
224
225
226 const auto geo_id_sl_map = makeGeoIdSourceLinkMap(tg, measurements);
227
228 auto hits = std::chrono::high_resolution_clock::now();
229 profiling_map_["hits"] +=
230 std::chrono::duration<double, std::milli>(hits - setup).count();
231
232
233
234
235 const std::vector<ldmx::Track> seed_tracks =
236 event.getCollection<
ldmx::Track>(seed_coll_name_, input_pass_name_);
237
238 ldmx_log(info) << "Number of " << seed_coll_name_
239 << " seed tracks = " << seed_tracks.size();
240
241 if (seed_tracks.empty()) {
242 std::vector<ldmx::Track> empty;
243 ldmx_log(info) << "No seed tracks, returning...";
244 event.add(out_trk_collection_, empty);
245 return;
246 }
247
248
249 std::vector<Acts::BoundTrackParameters> start_parameters;
250
251 ldmx_log(debug) << "Transform the seed track to bound parameters";
252 int seed_track_index{0};
253 for (auto& seed : seed_tracks) {
254
255 std::shared_ptr<Acts::PerigeeSurface> perigee_surface =
256 Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(
257 seed.getPerigeeX(), seed.getPerigeeY(), seed.getPerigeeZ()));
258
259 Acts::BoundVector param_vec;
260 param_vec << seed.getD0(), seed.getZ0(), seed.getPhi(), seed.getTheta(),
261 seed.getQoP(), seed.getT();
262
263 Acts::BoundSquareMatrix cov_mat =
264 tracking::sim::utils::unpackCov(seed.getPerigeeCov());
265
266 ldmx_log(debug) << " For seed index_ = " << seed_track_index
267 << ": Perigee X / Y / Z = " << seed.getPerigeeX() << " / "
268 << seed.getPerigeeY() << " / " << seed.getPerigeeZ()
269 << ", D0 = " << param_vec[0] << ", Z0 = " << param_vec[1]
270 << ", Phi = " << param_vec[2]
271 << ", Theta = " << param_vec[3]
272 << ", QoP = " << param_vec[4]
273 << ", Time = " << param_vec[5];
274
275 ldmx_log(debug) << " Cov matrix diagonal (" << cov_mat(0, 0) << ", "
276 << cov_mat(1, 1) << ", " << cov_mat(2, 2) << ")";
277
278
279 auto part_hypo{Acts::SinglyChargedParticleHypothesis::electron()};
280 start_parameters.push_back(Acts::BoundTrackParameters(
281 perigee_surface, param_vec, cov_mat, part_hypo));
282
283
285
286 seed_track_index++;
287 }
288
289 auto seeds = std::chrono::high_resolution_clock::now();
290 profiling_map_["seeds"] +=
291 std::chrono::duration<double, std::milli>(seeds - hits).count();
292
293 Acts::GainMatrixUpdater kf_updater;
294
295
296
297
298 Acts::MeasurementSelector::Config measurement_selector_cfg = {
299
300 {Acts::GeometryIdentifier(), {{}, {outlier_pval_}, {1u}}},
301 };
302
303 Acts::MeasurementSelector meas_sel{measurement_selector_cfg};
304
306
307 Acts::CombinatorialKalmanFilterExtensions<TrackContainer> ckf_extensions;
308
309 if (use1_dmeasurements_) {
310 ckf_extensions.calibrator
312 Acts::VectorMultiTrajectory>>(&calibrator);
313 } else {
314 ckf_extensions.calibrator
316 Acts::VectorMultiTrajectory>>(&calibrator);
317 }
318
319 ckf_extensions.updater.connect<
320 &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(
321 &kf_updater);
322
323 ckf_extensions.measurementSelector
324 .connect<&Acts::MeasurementSelector::select<Acts::VectorMultiTrajectory>>(
325 &meas_sel);
326
327 ldmx_log(debug) << "SourceLinkAccessor...";
328
329
330 struct SourceLinkAccIt {
331 using BaseIt = decltype(geo_id_sl_map.begin());
332 BaseIt it_;
333
334#pragma GCC diagnostic push
335#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
336
337 using difference_type = typename BaseIt::difference_type;
338 using iterator_category = typename BaseIt::iterator_category;
339
340 using value_type = Acts::SourceLink;
341 using pointer = typename BaseIt::pointer;
342 using reference = value_type&;
343#pragma GCC diagnostic pop
344
345 SourceLinkAccIt& operator++() {
346 ++it_;
347 return *this;
348 }
349 bool operator==(const SourceLinkAccIt& other) const {
350 return it_ == other.it_;
351 }
352 bool operator!=(const SourceLinkAccIt& other) const {
353 return !(*this == other);
354 }
355
356
357
358 value_type operator*() const { return value_type{it_->second}; }
359 };
360
361 auto source_link_accessor = [&](const Acts::Surface& surface)
362 -> std::pair<SourceLinkAccIt, SourceLinkAccIt> {
363 auto [begin, end] = geo_id_sl_map.equal_range(surface.geometryId());
364 return {SourceLinkAccIt{begin}, SourceLinkAccIt{end}};
365 };
366
367 Acts::SourceLinkAccessorDelegate<SourceLinkAccIt>
368 source_link_accessor_delegate;
369 source_link_accessor_delegate
370 .connect<&decltype(source_link_accessor)::operator(),
371 decltype(source_link_accessor)>(&source_link_accessor);
372
373 ldmx_log(debug) << "Setting up surfaces...";
374
375 std::shared_ptr<const Acts::PerigeeSurface> origin_surface =
376 Acts::Surface::makeShared<Acts::PerigeeSurface>(
377 Acts::Vector3(0., 0., 0.));
378
379 ldmx_log(debug) << "About to run CKF...";
380
381
382 auto ckf_setup = std::chrono::high_resolution_clock::now();
383 profiling_map_["ckf_setup"] +=
384 std::chrono::duration<double, std::milli>(ckf_setup - seeds).count();
385
386 Acts::VectorTrackContainer vtc;
387 Acts::VectorMultiTrajectory mtj;
388 Acts::TrackContainer tc{vtc, mtj};
389
390
391
392 ldmx_log(debug) << "Loop on the track candidates";
393 for (size_t track_id = 0u; track_id < start_parameters.size(); ++track_id) {
394 ldmx_log(debug) << "---------------------------";
395 ldmx_log(debug) << "Candidate Track ID = " << track_id;
396
397 const Acts::CombinatorialKalmanFilterOptions<SourceLinkAccIt,
398 TrackContainer>
399 ckf_options(TrackingGeometryUser::geometryContext(),
400 TrackingGeometryUser::magneticFieldContext(),
401 TrackingGeometryUser::calibrationContext(),
402 source_link_accessor_delegate, ckf_extensions,
403 propagator_options, true ,
404 false );
405
406 ldmx_log(debug) << " Checking options: multiple scattering = "
407 << ckf_options.multipleScattering
408 << " energy loss = " << ckf_options.energyLoss;
409 auto results =
410 ckf_->findTracks(start_parameters.at(track_id), ckf_options, tc);
411
412 auto start_params = start_parameters.at(track_id).parameters().transpose();
413 ldmx_log(debug)
414 << " Checking CKF success for track candidate with params: "
415 << " D0 = " << start_params[0] << " Z0 = " << start_params[1]
416 << ", Phi = " << start_params[2] << " Theta = " << start_params[3]
417 << ", QoP = " << start_params[4] << " Time = " << start_params[5];
418 if (not results.ok()) {
419 ldmx_log(debug) << " CKF failed!";
420 continue;
421 } else {
422 ldmx_log(debug) << " CKF succeded!";
423 }
424
425 auto& tracks_from_seed = results.value();
426 if (tracks_from_seed.size() != 1) {
427 ldmx_log(info) << " tracksFromSeed.size = " << tracks_from_seed.size();
428 }
429
430 for (auto& track : tracks_from_seed) {
431
432 Acts::smoothTrack(geometryContext(), track);
433
436
437
438 bool success = trk_extrap_->trackStateAtSurface(
439 track, target_surface_, ts_at_target, ldmx::TrackStateType::AtTarget);
440
441 if (success) {
442 ldmx_log(debug) << " Successfully obtained TrackState at target";
443
444 ldmx_log(debug) << " Parameters At Target: Loc0 = "
445 << ts_at_target.params_[0] << ", Loc1 "
446 << ts_at_target.params_[1]
447 << ", phi = " << ts_at_target.params_[2]
448 << ", theta = " << ts_at_target.params_[3] << ", QoP "
449 << ts_at_target.params_[4];
450
451 trk.addTrackState(ts_at_target);
452 } else {
453 ldmx_log(info) << " Could not extrapolate to target! nhits = "
454 << track.nMeasurements() << " Printing track states: ";
455 for (const auto ts : track.trackStatesReversed()) {
456 if (ts.hasSmoothed()) {
457 ldmx_log(info) << " Track momentum for smoothed track state = "
458 << 1 / ts.smoothed()[Acts::eBoundQOverP];
459 ldmx_log(info) << " Parameters: " << ts.smoothed().transpose();
460 } else {
461 ldmx_log(info) << " Track state not smoothed!";
462 }
463 }
464 ldmx_log(info) << " ...skipping this track candidate...";
465 continue;
466 }
467
468
469
470
471
472
473
474
475 Acts::BoundTrackParameters bound_state_at_target =
476 tracking::sim::utils::btp(ts_at_target, target_surface_, 11);
477 track.setReferenceSurface(target_surface_);
478 track.parameters() = bound_state_at_target.parameters();
479
480
481
482 const Acts::BoundVector& track_pars = track.parameters();
483
484
485 ldmx_log(debug) << " Track parameters (at target): Loc0 = "
486 << track_pars[Acts::eBoundLoc0]
487 << ", Loc1 = " << track_pars[Acts::eBoundLoc1]
488 << ", Theta = " << track_pars[Acts::eBoundTheta]
489 << ", Phi = " << track_pars[Acts::eBoundPhi]
490 << ", chi2 = " << track.chi2()
491 << ", nHits = " << track.nMeasurements();
492
493
494 trk.setPerigeeLocation(0, 0, 0);
496 trk.setPerigeeCov(ts_at_target.cov_);
497
498 trk.setChi2(track.chi2());
499 trk.setNhits(track.nMeasurements());
500
501
502 trk.setNdf(track.nMeasurements() - 5);
503 trk.setNsharedHits(track.nSharedHits());
504
505 ldmx_log(debug) << " Track momentum: px = " << track.momentum()[0]
506 << " py = " << track.momentum()[1]
507 << " pz = " << track.momentum()[2];
508 trk.setMomentum(track.momentum()[0], track.momentum()[1],
509 track.momentum()[2]);
510
511
512 if ((trk.getNhits() <= min_hits_) || (abs(1. / trk.getQoP()) <= 0.05)) {
513 ldmx_log(debug)
514 << " > Track candidate did NOT meet the requirements: Nhits = "
515 << trk.getNhits() << " and p = " << abs(1. / trk.getQoP())
516 << " GeV";
517
518
519 continue;
520 }
521
522
523 ldmx_log(debug) << " Add measurements to the final track from "
524 << track.nTrackStates() << " TrackStates with "
525 << track.nMeasurements() << " measurements";
526
527 int trk_state_index{0};
528 for (const auto ts : track.trackStatesReversed()) {
529
530 ldmx_log(debug) << " Checking Track State index_ = "
531 << trk_state_index << " at location "
532 << ts.referenceSurface()
533 .transform(geometryContext())
534 .translation()
535 .transpose();
536
537 if (ts.hasSmoothed()) {
538 ldmx_log(debug) << " Smoothed track parameters: "
539 << ts.smoothed().transpose();
540
541
542 }
543
544
545 auto type_flags = ts.typeFlags();
546
547 if (type_flags.test(Acts::TrackStateFlag::MeasurementFlag) &&
548 ts.hasUncalibratedSourceLink()) {
550 ts.getUncalibratedSourceLink()
551 .template get<acts_examples::IndexSourceLink>();
552
554 ldmx_log(debug) << " Adding measurement to ldmx::track with "
555 "source link index_ = "
557 ldmx_log(trace) << " Measurement:\n" << ldmx_meas;
558 trk.addMeasurementIndex(sl.
index());
559 } else {
560 ldmx_log(debug) << " This TrackState is not a measurement";
561 }
562 trk_state_index++;
563 }
564
565 ldmx_log(debug) << " Starting extrapolations";
566
567
568 const double ecal_scoring_plane = 240.5;
569 Acts::Vector3 pos(ecal_scoring_plane, 0., 0.);
570 Acts::Translation3 surf_translation(pos);
571 Acts::Transform3 surf_transform(surf_translation * surf_rotation_);
572 const std::shared_ptr<Acts::PlaneSurface> ecal_surface =
573 Acts::Surface::makeShared<Acts::PlaneSurface>(surf_transform);
574
575
576 const std::shared_ptr<Acts::Surface> beam_origin_surface =
577 tracking::sim::utils::unboundSurface(-700);
578
579 if (tagger_tracking_) {
580 ldmx_log(debug) << " Beam Origin Extrapolation";
582 success = trk_extrap_->trackStateAtSurface(
583 track, beam_origin_surface, ts_at_beam_origin,
584 ldmx::TrackStateType::AtBeamOrigin);
585
586 if (success) {
587 trk.addTrackState(ts_at_beam_origin);
588 ldmx_log(debug)
589 << " Successfully obtained TrackState at beam origin";
590 }
591 }
592
593
594 if (!tagger_tracking_) {
595 ldmx_log(debug) << " Ecal Extrapolation";
597 success = trk_extrap_->trackStateAtSurface(
598 track, ecal_surface, ts_at_ecal, ldmx::TrackStateType::AtECAL);
599
600 if (success) {
601 trk.addTrackState(ts_at_ecal);
602 ldmx_log(debug) << " Successfully obtained TrackState at Ecal";
603 ldmx_log(debug) << " Parameters At Ecal: Loc0 = "
604 << ts_at_ecal.params_[0]
605 << ", Loc1 = " << ts_at_ecal.params_[1]
606 << ", phi = " << ts_at_ecal.params_[2]
607 << ", theta = " << ts_at_ecal.params_[3]
608 << ", QoP = " << ts_at_ecal.params_[4];
609 } else {
610 ldmx_log(info) << " Could not extrapolate to ECAL!! Please check "
611 "the track states";
612 }
613 }
614
615
616 if (truth_matching_tool) {
617 auto truth_info = truth_matching_tool->truthMatch(trk);
618 trk.setTrackID(truth_info.track_id_);
619 trk.setPdgID(truth_info.pdg_id_);
620 trk.setTruthProb(truth_info.truth_prob_);
621 }
622
623
624 ldmx_log(debug)
625 << " > Adding the track candidate to the track collection";
626 tracks.push_back(trk);
627 ntracks_++;
628 }
629 }
630
631 ldmx_log(info) << "Number of CKF tracks " << tracks.size();
632
633 auto ckf_run = std::chrono::high_resolution_clock::now();
634 profiling_map_["ckf_run"] +=
635 std::chrono::duration<double, std::milli>(ckf_run - ckf_setup).count();
636
637
638 auto shared_hits = computeSharedHits(
639 tracks, measurements, tg, tracking::sim::utils::sourceLinkHash,
640 tracking::sim::utils::sourceLinkEquality);
641 for (std::size_t i_track = 0; i_track < shared_hits.size(); ++i_track) {
642 tracks[i_track].setNsharedHits(shared_hits[i_track].size());
643 for (auto idx : shared_hits[i_track]) {
644 tracks[i_track].addSharedIndex(idx);
645 }
646 }
647
648 auto result_loop = std::chrono::high_resolution_clock::now();
649 profiling_map_["result_loop"] +=
650 std::chrono::duration<double, std::milli>(result_loop - ckf_run).count();
651
652
653 event.add(out_trk_collection_, tracks);
654
655 auto end = std::chrono::high_resolution_clock::now();
656
657
658 auto diff = end - start;
659 processing_time_ += std::chrono::duration<double, std::milli>(diff).count();
660}
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 setPerigeeParameters(const std::vector< double > &par)
d_0 z_0 phi_0 theta q/p t
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.