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 truthMatchingTool_ = std::make_shared<tracking::sim::TruthMatchingTool>();
48 "out_seed_collection",
getName() +
"SeedTracks");
52 "input_hits_collection",
"TaggerSimHits");
56 "tagger_trks_collection",
"TaggerTracks");
59 "perigee_location", {-700, 0., 0.});
61 parameters.getParameter<
double>(
"pmin", 0.05 * Acts::UnitConstants::GeV);
62 pmax_ = parameters.getParameter<
double>(
"pmax", 8 * Acts::UnitConstants::GeV);
64 parameters.getParameter<
double>(
"d0max", -15. * Acts::UnitConstants::mm);
66 parameters.getParameter<
double>(
"d0min", -45. * Acts::UnitConstants::mm);
68 parameters.getParameter<
double>(
"z0max", 60. * Acts::UnitConstants::mm);
70 phicut_ = parameters.getParameter<
double>(
"phicut", 0.1);
71 thetacut_ = parameters.getParameter<
double>(
"thetacut", 0.2);
73 loc0cut_ = parameters.getParameter<
double>(
"loc0cut", 0.1);
74 loc1cut_ = parameters.getParameter<
double>(
"loc1cut", 0.3);
76 strategies_ = parameters.getParameter<std::vector<std::string>>(
77 "strategies", {
"0,1,2,3,4"});
79 inflate_factors_ = parameters.getParameter<std::vector<double>>(
80 "inflate_factors", {10., 10., 10., 10., 10., 10.});
82 bfield_ = parameters.getParameter<
double>(
"bfield", 1.5);
85 parameters.getParameter<std::string>(
"input_pass_name",
"");
91 auto start = std::chrono::high_resolution_clock::now();
92 ldmx::Tracks seed_tracks;
97 std::map<int, ldmx::SimParticle> particleMap;
99 const std::vector<ldmx::Measurement> measurements =
103 std::vector<ldmx::Track> tagger_tracks;
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 trackState = ts.value();
127 Acts::BoundSquareMatrix cov =
128 tracking::sim::utils::unpackCov(trackState.cov);
129 double locu = trackState.params[0];
130 double locv = trackState.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(geometry_context(), local_pos, dummy);
148 target_pseudo_meas.push_back(pseudo_meas);
152 if (event.
exists(
"SimParticles")) {
154 truthMatchingTool_->setup(particleMap, measurements);
157 ldmx_log(debug) <<
"Preparing the strategies";
164 std::vector<int> strategy = {0, 1, 2, 3, 4};
165 bool success = GroupStrips(measurements, strategy);
166 if (success) FindSeedsFromMap(seed_tracks, target_pseudo_meas);
180 ntracks_ += seed_tracks.size();
183 auto end = std::chrono::high_resolution_clock::now();
188 auto diff = end - start;
189 processing_time_ += std::chrono::duration<double, std::milli>(diff).count();
225 const ldmx::Measurements& vmeas,
double xOrigin,
226 const Acts::Vector3& perigee_location,
227 const ldmx::Measurements& pmeas_tgt) {
236 Acts::ActsMatrix<5, 5> A = Acts::ActsMatrix<5, 5>::Zero();
237 Acts::ActsVector<5> Y = Acts::ActsVector<5>::Zero();
239 for (
auto meas : vmeas) {
240 double xmeas = meas.getGlobalPosition()[0] - xOrigin;
243 const Acts::Surface* hit_surface = geometry().getSurface(meas.getLayerID());
246 auto rot = hit_surface->transform(geometry_context()).rotation();
247 auto tr = hit_surface->transform(geometry_context()).translation();
249 auto rotl2g = rot.transpose();
252 Acts::Vector2 loc{meas.getLocalPosition()[0], 0.};
254 xhit_.push_back(xmeas);
255 yhit_.push_back(meas.getGlobalPosition()[1]);
256 zhit_.push_back(meas.getGlobalPosition()[2]);
258 Acts::ActsMatrix<2, 5> A_i;
260 A_i(0, 0) = rotl2g(0, 1);
261 A_i(0, 1) = rotl2g(0, 1) * xmeas;
262 A_i(0, 2) = rotl2g(0, 1) * xmeas * xmeas;
263 A_i(0, 3) = rotl2g(0, 2);
264 A_i(0, 4) = rotl2g(0, 2) * xmeas;
266 A_i(1, 0) = rotl2g(1, 1);
267 A_i(1, 1) = rotl2g(1, 1) * xmeas;
268 A_i(1, 2) = rotl2g(1, 1) * xmeas * xmeas;
269 A_i(1, 3) = rotl2g(1, 2);
270 A_i(1, 4) = rotl2g(1, 2) * xmeas;
273 Acts::Vector2 offset = (rot.transpose() * tr).topRows<2>();
274 Acts::Vector2 xoffset = {rotl2g(0, 0) * xmeas, rotl2g(1, 0) * xmeas};
276 loc(0) = meas.getLocalPosition()[0];
278 double uError = sqrt(vmeas[0].getLocalCovariance()[0]);
281 double vError = 40. / sqrt(12);
283 Acts::ActsMatrix<2, 2> W_i =
284 Acts::ActsMatrix<2, 2>::Zero();
286 W_i(0, 0) = 1. / (uError * uError);
287 W_i(1, 1) = 1. / (vError * vError);
289 Acts::Vector2 Yprime_i = loc + offset - xoffset;
290 Y += (A_i.transpose()) * W_i * Yprime_i;
292 Acts::ActsMatrix<2, 5> WA_i = (W_i * A_i);
293 A += A_i.transpose() * WA_i;
296 Acts::ActsVector<5> B;
306 Acts::ActsVector<3> ref{0., 0., 0.};
308 double relativePerigeeX = perigee_location(0) - xOrigin;
310 std::shared_ptr<const Acts::PerigeeSurface> seed_perigee =
311 Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(
312 relativePerigeeX, perigee_location(1), perigee_location(2)));
315 Acts::Vector3 seed_pos{relativePerigeeX,
316 B(0) + B(1) * relativePerigeeX +
317 B(2) * relativePerigeeX * relativePerigeeX,
318 B(3) + B(4) * relativePerigeeX};
319 Acts::Vector3 dir{1, B(1) + 2 * B(2) * relativePerigeeX, B(4)};
324 0.3 * bfield_ * (1. / (2. * abs(B(2)))) * 0.001;
328 Acts::Vector3 seed_mom = p * dir / Acts::UnitConstants::MeV;
330 B(2) < 0 ? -1 * Acts::UnitConstants::e : +1 * Acts::UnitConstants::e;
346 (*seed_perigee).intersect(geometry_context(), seed_pos, dir);
348 Acts::FreeVector seed_free = tracking::sim::utils::toFreeParameters(
349 intersection.intersections()[0].position(), seed_mom, q);
351 auto bound_params = Acts::transformFreeToBoundParameters(
352 seed_free, *seed_perigee, geometry_context())
355 ldmx_log(trace) <<
"bound parameters at perigee location" << bound_params;
357 Acts::BoundVector stddev;
359 double sigma_p = 0.75 * p * Acts::UnitConstants::GeV;
360 stddev[Acts::eBoundLoc0] =
361 inflate_factors_[Acts::eBoundLoc0] * 2 * Acts::UnitConstants::mm;
362 stddev[Acts::eBoundLoc1] =
363 inflate_factors_[Acts::eBoundLoc1] * 5 * Acts::UnitConstants::mm;
364 stddev[Acts::eBoundPhi] =
365 inflate_factors_[Acts::eBoundPhi] * 5 * Acts::UnitConstants::degree;
366 stddev[Acts::eBoundTheta] =
367 inflate_factors_[Acts::eBoundTheta] * 5 * Acts::UnitConstants::degree;
368 stddev[Acts::eBoundQOverP] =
369 inflate_factors_[Acts::eBoundQOverP] * (1. / p) * (1. / p) * sigma_p;
370 stddev[Acts::eBoundTime] =
371 inflate_factors_[Acts::eBoundTime] * 1000 * Acts::UnitConstants::ns;
374 <<
"Making covariance matrix as diagonal matrix with inflated terms";
375 Acts::BoundSquareMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal();
377 ldmx_log(debug) <<
"...now putting together the seed track ...";
380 trk.setPerigeeLocation(perigee_location(0), perigee_location(1),
381 perigee_location(2));
385 trk.setNsharedHits(0);
386 std::vector<double> v_seed_params(
387 (bound_params).data(),
388 bound_params.data() + bound_params.rows() * bound_params.cols());
389 std::vector<double> v_seed_cov;
390 tracking::sim::utils::flatCov(bound_cov, v_seed_cov);
392 trk.setPerigeeCov(v_seed_cov);
397 trk.setPosition(seed_free[Acts::eFreePos0], seed_free[Acts::eFreePos1],
398 seed_free[Acts::eFreePos2]);
399 trk.setMomentum(seed_free[Acts::eFreeDir0], seed_free[Acts::eFreeDir1],
400 seed_free[Acts::eFreeDir2]);
403 <<
"...making the ParticleHypothesis ...assume electron for now";
404 auto partHypo{Acts::SinglyChargedParticleHypothesis::electron()};
406 ldmx_log(debug) <<
"Making BoundTrackParameters seedParameters";
407 Acts::BoundTrackParameters seedParameters(
408 seed_perigee, std::move(bound_params), bound_cov, partHypo);
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(ldmx::Tracks& 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();
484 while (it[0] != groups_iter->second.end()) {
496 std::vector<ldmx::Measurement> meas_for_seeds;
497 meas_for_seeds.reserve(5);
499 ldmx_log(debug) <<
" Grouping ";
501 for (
int j = 0; j < K; j++) {
503 meas_for_seeds.push_back(*meas);
506 std::sort(meas_for_seeds.begin(), meas_for_seeds.end(),
508 return m1.getGlobalPosition()[0] < m2.getGlobalPosition()[0];
511 if (meas_for_seeds.size() < 5) {
516 ldmx_log(debug) <<
"making seedTrack";
522 SeedTracker(meas_for_seeds, meas_for_seeds.at(2).getGlobalPosition()[0],
528 if (1. / abs(seedTrack.getQoP()) <
pmin_) {
531 }
else if (1. / abs(seedTrack.getQoP()) >
pmax_) {
539 else if (abs(seedTrack.getZ0()) >
z0max_) {
542 }
else if (seedTrack.getD0() <
d0min_) {
545 }
else if (seedTrack.getD0() >
d0max_) {
548 }
else if (abs(seedTrack.getPhi()) >
phicut_) {
551 }
else if (abs(seedTrack.getTheta() - piover2_) >
thetacut_) {
562 if (pmeas.size() > 0) {
569 for (
auto tgt_pseudomeas : pmeas) {
573 seedTrack.getD0() - tgt_pseudomeas.getLocalPosition()[0];
575 seedTrack.getZ0() - tgt_pseudomeas.getLocalPosition()[1];
577 if (abs(delta_loc0) <
loc0cut_ && abs(delta_loc1) < loc1cut_) {
586 if (truthMatchingTool_->configured()) {
587 auto truthInfo = truthMatchingTool_->TruthMatch(meas_for_seeds);
588 seedTrack.setTrackID(truthInfo.trackID);
589 seedTrack.setPdgID(truthInfo.pdgID);
590 seedTrack.setTruthProb(truthInfo.truthProb);
593 seeds.push_back(seedTrack);
605 ldmx_log(debug) <<
"Go to the next combination";
609 (i > 0) && (it[i] == (std::next(groups_iter, i))->second.end()); --i) {
610 it[i] = std::next(groups_iter, i)->second.begin();
#define DECLARE_PRODUCER_NS(NS, 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.
void setTime(const float &t)
Set the measurement time in ns.
void setGlobalPosition(const float &x, const float &y, const float &z)
Set the global position i.e.
void setLocalPosition(const float &u, const float &v)
Set the local position i.e.
void setLocalCovariance(const float &cov_uu, const float &cov_vv)
Set cov(U,U) and cov(V, V).
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...