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 std::vector<ldmx::Track> seed_tracks;
96 std::map<int, ldmx::SimParticle> particle_map;
101 std::vector<ldmx::Track> tagger_tracks;
103 tagger_trks_event_collection_passname_)) {
109 std::shared_ptr<Acts::Surface> tgt_surf =
110 tracking::sim::utils::unboundSurface(0.);
114 ldmx::Measurements target_pseudo_meas;
116 for (
auto tagtrk : tagger_tracks) {
126 const auto& perigee_cov = tagtrk.getPerigeeCov();
127 if (!perigee_cov.empty()) {
128 Acts::BoundMatrix cov = tracking::sim::utils::unpackCov(perigee_cov);
129 double locu = tagtrk.getD0();
130 double locv = tagtrk.getZ0();
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::Matrix<5, 5> a = Acts::Matrix<5, 5>::Zero();
239 Acts::Vector<5> y = Acts::Vector<5>::Zero();
241 for (
auto meas : vmeas) {
242 double xmeas = meas.getGlobalPosition()[0] - xOrigin;
245 const Acts::Surface* hit_surface = geometry().getSurface(meas.getLayerID());
249 hit_surface->localToGlobalTransform(geometryContext()).rotation();
251 hit_surface->localToGlobalTransform(geometryContext()).translation();
253 auto rotl2g = rot.transpose();
256 Acts::Vector2 loc{meas.getLocalPosition()[0], 0.};
258 xhit_.push_back(xmeas);
259 yhit_.push_back(meas.getGlobalPosition()[1]);
260 zhit_.push_back(meas.getGlobalPosition()[2]);
262 Acts::Matrix<2, 5> a_i;
264 a_i(0, 0) = rotl2g(0, 1);
265 a_i(0, 1) = rotl2g(0, 1) * xmeas;
266 a_i(0, 2) = rotl2g(0, 1) * xmeas * xmeas;
267 a_i(0, 3) = rotl2g(0, 2);
268 a_i(0, 4) = rotl2g(0, 2) * xmeas;
270 a_i(1, 0) = rotl2g(1, 1);
271 a_i(1, 1) = rotl2g(1, 1) * xmeas;
272 a_i(1, 2) = rotl2g(1, 1) * xmeas * xmeas;
273 a_i(1, 3) = rotl2g(1, 2);
274 a_i(1, 4) = rotl2g(1, 2) * xmeas;
277 Acts::Vector2 offset = (rot.transpose() * tr).topRows<2>();
278 Acts::Vector2 xoffset = {rotl2g(0, 0) * xmeas, rotl2g(1, 0) * xmeas};
280 loc(0) = meas.getLocalPosition()[0];
283 Acts::Matrix<2, 2> w_i = Acts::Matrix<2, 2>::Zero();
288 Acts::Vector2 yprime_i = loc + offset - xoffset;
289 y += (a_i.transpose()) * w_i * yprime_i;
291 Acts::Matrix<2, 5> wa_i = (w_i * a_i);
292 a += a_i.transpose() * wa_i;
305 Acts::Vector<3> ref{0., 0., 0.};
311 double relative_perigee_x = perigee_location(0) - xOrigin;
313 std::shared_ptr<const Acts::PerigeeSurface> seed_perigee =
314 Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(
315 perigee_location(0), perigee_location(1), perigee_location(2)));
319 Acts::Vector3 seed_pos{perigee_location(0),
320 b(0) + b(1) * relative_perigee_x +
321 b(2) * relative_perigee_x * relative_perigee_x,
322 b(3) + b(4) * relative_perigee_x};
323 Acts::Vector3 dir{1, b(1) + 2 * b(2) * relative_perigee_x, b(4)};
328 double p = 0.3 * bfield_ * (1. / (2. * abs(b(2)))) * 0.001;
332 Acts::Vector3 seed_mom = p * dir / Acts::UnitConstants::MeV;
334 b(2) < 0 ? -1 * Acts::UnitConstants::e : +1 * Acts::UnitConstants::e;
350 (*seed_perigee).intersect(geometryContext(), seed_pos, dir);
352 Acts::FreeVector seed_free = tracking::sim::utils::toFreeParameters(
353 intersection[0].position(), seed_mom, q);
355 auto bound_params = Acts::transformFreeToBoundParameters(
356 seed_free, *seed_perigee, geometryContext())
359 ldmx_log(trace) <<
"bound parameters at perigee location" << bound_params;
361 Acts::BoundVector stddev;
363 double sigma_p = 0.75 * p * Acts::UnitConstants::GeV;
364 stddev[Acts::eBoundLoc0] =
365 inflate_factors_[Acts::eBoundLoc0] * 2 * Acts::UnitConstants::mm;
366 stddev[Acts::eBoundLoc1] =
367 inflate_factors_[Acts::eBoundLoc1] * 5 * Acts::UnitConstants::mm;
368 stddev[Acts::eBoundPhi] =
369 inflate_factors_[Acts::eBoundPhi] * 5 * Acts::UnitConstants::degree;
370 stddev[Acts::eBoundTheta] =
371 inflate_factors_[Acts::eBoundTheta] * 5 * Acts::UnitConstants::degree;
372 stddev[Acts::eBoundQOverP] =
373 inflate_factors_[Acts::eBoundQOverP] * (1. / p) * (1. / p) * sigma_p;
374 stddev[Acts::eBoundTime] =
375 inflate_factors_[Acts::eBoundTime] * 1000 * Acts::UnitConstants::ns;
378 <<
"Making covariance matrix as diagonal matrix with inflated terms";
379 Acts::BoundMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal();
381 ldmx_log(debug) <<
"...now putting together the seed track ...";
386 Acts::Vector3 perigee_ldmx =
387 tracking::sim::utils::acts2Ldmx(perigee_location);
388 trk.setPerigeeLocation(perigee_ldmx(0), perigee_ldmx(1), perigee_ldmx(2));
392 trk.setNsharedHits(0);
393 trk.setCharge(q < 0 ? -1 : 1);
394 std::vector<double> v_seed_params(
395 (bound_params).data(),
396 bound_params.data() + bound_params.rows() * bound_params.cols());
397 std::vector<double> v_seed_cov;
398 tracking::sim::utils::flatCov(bound_cov, v_seed_cov);
399 trk.setPerigeeParameters(v_seed_params);
400 trk.setPerigeeCov(v_seed_cov);
403 <<
"...making the ParticleHypothesis ...assume electron for now";
404 auto part_hypo{Acts::ParticleHypothesis::electron()};
406 ldmx_log(debug) <<
"Making BoundTrackParameters seedParameters";
407 Acts::BoundTrackParameters seed_parameters(
408 seed_perigee, std::move(bound_params), bound_cov, part_hypo);
410 ldmx_log(debug) <<
"Returning seed track";
418 ldmx_log(info) <<
"AVG Time/Event: " << std::fixed << std::setprecision(1)
419 << processing_time_ / nevents_ <<
" ms";
420 ldmx_log(info) <<
"Total Seeds/Events: " << ntracks_ <<
"/" << nevents_;
421 ldmx_log(info) <<
"Seeds discarded due to multiple hits on layers "
423 ldmx_log(info) <<
"not enough seed points " << nmissing_;
424 ldmx_log(info) <<
" nfailpmin=" << nfailpmin_;
425 ldmx_log(info) <<
" nfailpmax=" << nfailpmax_;
426 ldmx_log(info) <<
" nfaild0max=" << nfaild0max_;
427 ldmx_log(info) <<
" nfaild0min=" << nfaild0min_;
428 ldmx_log(info) <<
" nfailphicut=" << nfailphi_;
429 ldmx_log(info) <<
" nfailthetacut=" << nfailtheta_;
430 ldmx_log(info) <<
" nfailz0max=" << nfailz0max_;
437bool SeedFinderProcessor::groupStrips(
438 const std::vector<ldmx::Measurement>& measurements,
439 const std::vector<int> strategy) {
446 for (
auto& meas : measurements) {
447 ldmx_log(trace) << meas;
449 if (std::find(strategy.begin(), strategy.end(), meas.getLayer()) !=
451 ldmx_log(debug) <<
"Adding measurement from layer_ = " << meas.getLayer();
452 groups_map_[meas.getLayer()].push_back(&meas);
457 if (groups_map_.size() < 5)
467void SeedFinderProcessor::findSeedsFromMap(std::vector<ldmx::Track>& seeds,
468 const ldmx::Measurements& pmeas) {
469 std::map<int, std::vector<const ldmx::Measurement*>>::iterator groups_iter =
472 constexpr size_t k = 5;
473 std::vector<std::vector<const ldmx::Measurement*>::iterator> it;
476 unsigned int ikey = 0;
477 for (
auto& key : groups_map_) {
478 it[ikey] = key.second.begin();
485 while (it[0] != groups_iter->second.end()) {
497 std::vector<ldmx::Measurement> meas_for_seeds;
498 meas_for_seeds.reserve(5);
500 ldmx_log(debug) <<
" Grouping ";
502 for (
int j = 0; j < k; j++) {
504 meas_for_seeds.push_back(*meas);
507 std::sort(meas_for_seeds.begin(), meas_for_seeds.end(),
509 return m1.getGlobalPosition()[0] < m2.getGlobalPosition()[0];
512 if (meas_for_seeds.size() < 5) {
517 ldmx_log(debug) <<
"making seedTrack";
523 seedTracker(meas_for_seeds, meas_for_seeds.at(2).getGlobalPosition()[0],
529 if (1. / abs(seed_track.getQoP()) <
pmin_) {
532 }
else if (1. / abs(seed_track.getQoP()) >
pmax_) {
540 else if (abs(seed_track.getZ0()) >
z0max_) {
543 }
else if (seed_track.getD0() <
d0min_) {
546 }
else if (seed_track.getD0() >
d0max_) {
549 }
else if (abs(seed_track.getPhi()) >
phicut_) {
552 }
else if (abs(seed_track.getTheta() - piover2_) >
thetacut_) {
563 if (pmeas.size() > 0) {
570 for (
auto tgt_pseudomeas : pmeas) {
574 seed_track.getD0() - tgt_pseudomeas.getLocalPosition()[0];
576 seed_track.getZ0() - tgt_pseudomeas.getLocalPosition()[1];
578 if (abs(delta_loc0) <
loc0cut_ && abs(delta_loc1) < loc1cut_) {
587 if (truth_matching_tool_->configured()) {
588 auto truth_info = truth_matching_tool_->truthMatch(meas_for_seeds);
589 seed_track.setTrackID(truth_info.track_id_);
590 seed_track.setPdgID(truth_info.pdg_id_);
591 seed_track.setTruthProb(truth_info.truth_prob_);
594 seeds.push_back(seed_track);
606 ldmx_log(debug) <<
"Go to the next combination";
610 (i > 0) && (it[i] == (std::next(groups_iter, i))->second.end()); --i) {
611 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.
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...