1#include "Tracking/Reco/SeedFinderProcessor.h"
3#include "Acts/Definitions/TrackParametrization.hpp"
4#include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
6#include "Tracking/Sim/TrackingUtils.h"
42 truth_matching_tool_ = std::make_shared<tracking::sim::TruthMatchingTool>();
52 parameters.
get<std::string>(
"input_hits_collection",
"TaggerSimHits");
56 parameters.
get<std::string>(
"tagger_trks_collection",
"TaggerTracks");
59 parameters.
get<std::vector<double>>(
"perigee_location", {-700, 0., 0.});
60 pmin_ = parameters.
get<
double>(
"pmin", 0.05 * Acts::UnitConstants::GeV);
61 pmax_ = parameters.
get<
double>(
"pmax", 8 * Acts::UnitConstants::GeV);
62 d0max_ = parameters.
get<
double>(
"d0max", -15. * Acts::UnitConstants::mm);
63 d0min_ = parameters.
get<
double>(
"d0min", -45. * Acts::UnitConstants::mm);
64 z0max_ = parameters.
get<
double>(
"z0max", 60. * Acts::UnitConstants::mm);
65 phicut_ = parameters.
get<
double>(
"phicut", 0.1);
68 loc1cut_ = parameters.
get<
double>(
"loc1cut", 0.3);
70 parameters.
get<std::vector<std::string>>(
"strategies", {
"0,1,2,3,4"});
71 inflate_factors_ = parameters.
get<std::vector<double>>(
72 "inflate_factors", {10., 10., 10., 10., 10., 10.});
73 bfield_ = parameters.
get<
double>(
"bfield", 1.5);
74 input_pass_name_ = parameters.
get<std::string>(
"input_pass_name");
75 sim_particles_coll_name_ =
76 parameters.
get<std::string>(
"sim_particles_coll_name");
77 sim_particles_passname_ =
78 parameters.
get<std::string>(
"sim_particles_passname");
79 tagger_trks_event_collection_passname_ =
80 parameters.
get<std::string>(
"tagger_trks_event_collection_passname");
81 sim_particles_event_passname_ =
82 parameters.
get<std::string>(
"sim_particles_event_passname");
90 auto start = std::chrono::high_resolution_clock::now();
91 ldmx::Tracks seed_tracks;
96 std::map<int, ldmx::SimParticle> particle_map;
98 const std::vector<ldmx::Measurement> measurements =
102 std::vector<ldmx::Track> tagger_tracks;
104 tagger_trks_event_collection_passname_)) {
110 std::shared_ptr<Acts::Surface> tgt_surf =
111 tracking::sim::utils::unboundSurface(0.);
115 ldmx::Measurements target_pseudo_meas;
117 for (
auto tagtrk : tagger_tracks) {
118 auto ts = tagtrk.getTrackState(ldmx::TrackStateType::AtTarget);
124 if (ts.has_value()) {
125 auto track_state = ts.value();
127 Acts::BoundSquareMatrix cov =
128 tracking::sim::utils::unpackCov(track_state.cov_);
129 double locu = track_state.params_[0];
130 double locv = track_state.params_[1];
132 cov(Acts::BoundIndices::eBoundLoc0, Acts::BoundIndices::eBoundLoc0);
134 cov(Acts::BoundIndices::eBoundLoc1, Acts::BoundIndices::eBoundLoc1);
138 Acts::Vector3 dummy{0., 0., 0.};
139 Acts::Vector2 local_pos{locu, locv};
140 Acts::Vector3 global_pos =
141 tgt_surf->localToGlobal(geometryContext(), local_pos, dummy);
148 target_pseudo_meas.push_back(pseudo_meas);
152 if (event.
exists(sim_particles_coll_name_, sim_particles_event_passname_)) {
154 sim_particles_coll_name_, sim_particles_passname_);
155 truth_matching_tool_->setup(particle_map, measurements);
158 ldmx_log(debug) <<
"Preparing the strategies";
165 std::vector<int> strategy = {0, 1, 2, 3, 4};
166 bool success = groupStrips(measurements, strategy);
167 if (success) findSeedsFromMap(seed_tracks, target_pseudo_meas);
181 ntracks_ += seed_tracks.size();
184 auto end = std::chrono::high_resolution_clock::now();
189 auto diff = end - start;
190 processing_time_ += std::chrono::duration<double, std::milli>(diff).count();
227 const ldmx::Measurements& vmeas,
double xOrigin,
228 const Acts::Vector3& perigee_location,
229 const ldmx::Measurements& pmeas_tgt) {
238 Acts::ActsMatrix<5, 5> a = Acts::ActsMatrix<5, 5>::Zero();
239 Acts::ActsVector<5> y = Acts::ActsVector<5>::Zero();
241 for (
auto meas : vmeas) {
242 double xmeas = meas.getGlobalPosition()[0] - xOrigin;
245 const Acts::Surface* hit_surface = geometry().getSurface(meas.getLayerID());
248 auto rot = hit_surface->transform(geometryContext()).rotation();
249 auto tr = hit_surface->transform(geometryContext()).translation();
251 auto rotl2g = rot.transpose();
254 Acts::Vector2 loc{meas.getLocalPosition()[0], 0.};
256 xhit_.push_back(xmeas);
257 yhit_.push_back(meas.getGlobalPosition()[1]);
258 zhit_.push_back(meas.getGlobalPosition()[2]);
260 Acts::ActsMatrix<2, 5> a_i;
262 a_i(0, 0) = rotl2g(0, 1);
263 a_i(0, 1) = rotl2g(0, 1) * xmeas;
264 a_i(0, 2) = rotl2g(0, 1) * xmeas * xmeas;
265 a_i(0, 3) = rotl2g(0, 2);
266 a_i(0, 4) = rotl2g(0, 2) * xmeas;
268 a_i(1, 0) = rotl2g(1, 1);
269 a_i(1, 1) = rotl2g(1, 1) * xmeas;
270 a_i(1, 2) = rotl2g(1, 1) * xmeas * xmeas;
271 a_i(1, 3) = rotl2g(1, 2);
272 a_i(1, 4) = rotl2g(1, 2) * xmeas;
275 Acts::Vector2 offset = (rot.transpose() * tr).topRows<2>();
276 Acts::Vector2 xoffset = {rotl2g(0, 0) * xmeas, rotl2g(1, 0) * xmeas};
278 loc(0) = meas.getLocalPosition()[0];
281 Acts::ActsMatrix<2, 2> w_i = Acts::ActsMatrix<2, 2>::Zero();
286 Acts::Vector2 yprime_i = loc + offset - xoffset;
287 y += (a_i.transpose()) * w_i * yprime_i;
289 Acts::ActsMatrix<2, 5> wa_i = (w_i * a_i);
290 a += a_i.transpose() * wa_i;
293 Acts::ActsVector<5> b;
303 Acts::ActsVector<3> ref{0., 0., 0.};
305 double relative_perigee_x = perigee_location(0) - xOrigin;
307 std::shared_ptr<const Acts::PerigeeSurface> seed_perigee =
308 Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(
309 relative_perigee_x, perigee_location(1), perigee_location(2)));
312 Acts::Vector3 seed_pos{relative_perigee_x,
313 b(0) + b(1) * relative_perigee_x +
314 b(2) * relative_perigee_x * relative_perigee_x,
315 b(3) + b(4) * relative_perigee_x};
316 Acts::Vector3 dir{1, b(1) + 2 * b(2) * relative_perigee_x, b(4)};
321 double p = 0.3 * bfield_ * (1. / (2. * abs(b(2)))) * 0.001;
325 Acts::Vector3 seed_mom = p * dir / Acts::UnitConstants::MeV;
327 b(2) < 0 ? -1 * Acts::UnitConstants::e : +1 * Acts::UnitConstants::e;
343 (*seed_perigee).intersect(geometryContext(), seed_pos, dir);
345 Acts::FreeVector seed_free = tracking::sim::utils::toFreeParameters(
346 intersection.intersections()[0].position(), seed_mom, q);
348 auto bound_params = Acts::transformFreeToBoundParameters(
349 seed_free, *seed_perigee, geometryContext())
352 ldmx_log(trace) <<
"bound parameters at perigee location" << bound_params;
354 Acts::BoundVector stddev;
356 double sigma_p = 0.75 * p * Acts::UnitConstants::GeV;
357 stddev[Acts::eBoundLoc0] =
358 inflate_factors_[Acts::eBoundLoc0] * 2 * Acts::UnitConstants::mm;
359 stddev[Acts::eBoundLoc1] =
360 inflate_factors_[Acts::eBoundLoc1] * 5 * Acts::UnitConstants::mm;
361 stddev[Acts::eBoundPhi] =
362 inflate_factors_[Acts::eBoundPhi] * 5 * Acts::UnitConstants::degree;
363 stddev[Acts::eBoundTheta] =
364 inflate_factors_[Acts::eBoundTheta] * 5 * Acts::UnitConstants::degree;
365 stddev[Acts::eBoundQOverP] =
366 inflate_factors_[Acts::eBoundQOverP] * (1. / p) * (1. / p) * sigma_p;
367 stddev[Acts::eBoundTime] =
368 inflate_factors_[Acts::eBoundTime] * 1000 * Acts::UnitConstants::ns;
371 <<
"Making covariance matrix as diagonal matrix with inflated terms";
372 Acts::BoundSquareMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal();
374 ldmx_log(debug) <<
"...now putting together the seed track ...";
377 trk.setPerigeeLocation(perigee_location(0), perigee_location(1),
378 perigee_location(2));
382 trk.setNsharedHits(0);
383 std::vector<double> v_seed_params(
384 (bound_params).data(),
385 bound_params.data() + bound_params.rows() * bound_params.cols());
386 std::vector<double> v_seed_cov;
387 tracking::sim::utils::flatCov(bound_cov, v_seed_cov);
389 trk.setPerigeeCov(v_seed_cov);
394 trk.setPosition(seed_free[Acts::eFreePos0], seed_free[Acts::eFreePos1],
395 seed_free[Acts::eFreePos2]);
396 trk.setMomentum(seed_free[Acts::eFreeDir0], seed_free[Acts::eFreeDir1],
397 seed_free[Acts::eFreeDir2]);
400 <<
"...making the ParticleHypothesis ...assume electron for now";
401 auto part_hypo{Acts::SinglyChargedParticleHypothesis::electron()};
403 ldmx_log(debug) <<
"Making BoundTrackParameters seedParameters";
404 Acts::BoundTrackParameters seed_parameters(
405 seed_perigee, std::move(bound_params), bound_cov, part_hypo);
407 ldmx_log(debug) <<
"Returning seed track";
415 ldmx_log(info) <<
"AVG Time/Event: " << std::fixed << std::setprecision(1)
416 << processing_time_ / nevents_ <<
" ms";
417 ldmx_log(info) <<
"Total Seeds/Events: " << ntracks_ <<
"/" << nevents_;
418 ldmx_log(info) <<
"Seeds discarded due to multiple hits on layers "
420 ldmx_log(info) <<
"not enough seed points " << nmissing_;
421 ldmx_log(info) <<
" nfailpmin=" << nfailpmin_;
422 ldmx_log(info) <<
" nfailpmax=" << nfailpmax_;
423 ldmx_log(info) <<
" nfaild0max=" << nfaild0max_;
424 ldmx_log(info) <<
" nfaild0min=" << nfaild0min_;
425 ldmx_log(info) <<
" nfailphicut=" << nfailphi_;
426 ldmx_log(info) <<
" nfailthetacut=" << nfailtheta_;
427 ldmx_log(info) <<
" nfailz0max=" << nfailz0max_;
434bool SeedFinderProcessor::groupStrips(
435 const std::vector<ldmx::Measurement>& measurements,
436 const std::vector<int> strategy) {
443 for (
auto& meas : measurements) {
444 ldmx_log(trace) << meas;
446 if (std::find(strategy.begin(), strategy.end(), meas.getLayer()) !=
448 ldmx_log(debug) <<
"Adding measurement from layer_ = " << meas.getLayer();
449 groups_map_[meas.getLayer()].push_back(&meas);
454 if (groups_map_.size() < 5)
464void SeedFinderProcessor::findSeedsFromMap(ldmx::Tracks& seeds,
465 const ldmx::Measurements& pmeas) {
466 std::map<int, std::vector<const ldmx::Measurement*>>::iterator groups_iter =
469 constexpr size_t k = 5;
470 std::vector<std::vector<const ldmx::Measurement*>::iterator> it;
473 unsigned int ikey = 0;
474 for (
auto& key : groups_map_) {
475 it[ikey] = key.second.begin();
482 while (it[0] != groups_iter->second.end()) {
494 std::vector<ldmx::Measurement> meas_for_seeds;
495 meas_for_seeds.reserve(5);
497 ldmx_log(debug) <<
" Grouping ";
499 for (
int j = 0; j < k; j++) {
501 meas_for_seeds.push_back(*meas);
504 std::sort(meas_for_seeds.begin(), meas_for_seeds.end(),
506 return m1.getGlobalPosition()[0] < m2.getGlobalPosition()[0];
509 if (meas_for_seeds.size() < 5) {
514 ldmx_log(debug) <<
"making seedTrack";
520 seedTracker(meas_for_seeds, meas_for_seeds.at(2).getGlobalPosition()[0],
526 if (1. / abs(seed_track.getQoP()) <
pmin_) {
529 }
else if (1. / abs(seed_track.getQoP()) >
pmax_) {
537 else if (abs(seed_track.getZ0()) >
z0max_) {
540 }
else if (seed_track.getD0() <
d0min_) {
543 }
else if (seed_track.getD0() >
d0max_) {
546 }
else if (abs(seed_track.getPhi()) >
phicut_) {
549 }
else if (abs(seed_track.getTheta() - piover2_) >
thetacut_) {
560 if (pmeas.size() > 0) {
567 for (
auto tgt_pseudomeas : pmeas) {
571 seed_track.getD0() - tgt_pseudomeas.getLocalPosition()[0];
573 seed_track.getZ0() - tgt_pseudomeas.getLocalPosition()[1];
575 if (abs(delta_loc0) <
loc0cut_ && abs(delta_loc1) < loc1cut_) {
584 if (truth_matching_tool_->configured()) {
585 auto truth_info = truth_matching_tool_->truthMatch(meas_for_seeds);
586 seed_track.setTrackID(truth_info.track_id_);
587 seed_track.setPdgID(truth_info.pdg_id_);
588 seed_track.setTruthProb(truth_info.truth_prob_);
591 seeds.push_back(seed_track);
603 ldmx_log(debug) <<
"Go to the next combination";
607 (i > 0) && (it[i] == (std::next(groups_iter, i))->second.end()); --i) {
608 it[i] = std::next(groups_iter, i)->second.begin();
#define DECLARE_PRODUCER(CLASS)
Macro which allows the framework to construct a producer given its name during configuration.
std::string getName() const
Get the processor name.
Implements an event buffer system for storing event data.
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 which represents the process under execution.
Class encapsulating parameters for configuring a processor.
const T & get(const std::string &name) const
Retrieve the parameter of the given name.
void setLocalPosition(const float &meas_u, const float &meas_v)
Set the local position i.e.
void setGlobalPosition(const float &meas_x, const float &meas_y, const float &meas_z)
Set the global position i.e.
void setLocalCovariance(const float &cov_uu, const float &cov_vv)
Set cov(U,U) and cov(V, V).
void setTime(const float &meas_t)
Set the measurement time in ns.
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
double thetacut_
ThetaRange.
double loc0cut_
loc0 / loc1 cuts
SeedFinderProcessor(const std::string &name, framework::Process &process)
Constructor.
std::string out_seed_collection_
The name of the output collection of seeds to be stored.
double pmax_
Maximum cut on the momentum of the seeds.
std::vector< std::string > strategies_
List of stragies for seed finding.
std::string input_hits_collection_
The name of the input hits collection to use in finding seeds..
void onProcessEnd() override
Callback for the EventProcessor to take any necessary action when the processing of events finishes,...
void produce(framework::Event &event) override
Run the processor and create a collection of results which indicate if a charge particle can be found...
void onProcessStart() override
Callback for the EventProcessor to take any necessary action when the processing of events starts,...
double pmin_
Minimum cut on the momentum of the seeds.
std::string tagger_trks_collection_
The name of the tagger Tracks (only for Recoil Seeding)
std::vector< double > perigee_location_
Location of the perigee for the helix track parameters.
double d0max_
Max d0 allowed for the seeds.
void configure(framework::config::Parameters ¶meters) override
Configure the processor using the given user specified parameters.
double d0min_
Min d0 allowed for the seeds.
double z0max_
Max z0 allowed for the seeds.
a helper base class providing some methods to shorten access to common conditions used within the tra...
The measurement calibrator can be a function or a class/struct able to retrieve the sim hits containe...