Public Member Functions
	LinearSeedFinder (const std::string &name, framework::Process &process)
	Constructor.

virtual	~LinearSeedFinder ()=default
	Destructor.

void	onProcessStart () override
	Setup the truth matching.

void	onProcessEnd () override
	Output event statistics.

void	configure (framework::config::Parameters &parameters) override
	Configure the processor using the given user specified parameters.

void	produce (framework::Event &event) override
	Run the processor and create a collection of results which indicate if a charge particle can be found by the recoil tracker.

Public Member Functions inherited from tracking::reco::TrackingGeometryUser
	TrackingGeometryUser (const std::string &name, framework::Process &p)

Public Member Functions inherited from framework::Producer
	Producer (const std::string &name, Process &process)
	Class constructor.

virtual void	process (Event &event) final
	Processing an event for a Producer is calling produce.

Public Member Functions inherited from framework::EventProcessor
	DECLARE_FACTORY (EventProcessor, EventProcessor *, const std::string &, Process &)
	declare that we have a factory for this class

	EventProcessor (const std::string &name, Process &process)
	Class constructor.

virtual	~EventProcessor ()=default
	Class destructor.

virtual void	beforeNewRun (ldmx::RunHeader &run_header)
	Callback for Producers to add parameters to the run header before conditions are initialized.

virtual void	onNewRun (const ldmx::RunHeader &run_header)
	Callback for the EventProcessor to take any necessary action when the run being processed changes.

virtual void	onFileOpen (EventFile &event_file)
	Callback for the EventProcessor to take any necessary action when a new event input ROOT file is opened.

virtual void	onFileClose (EventFile &event_file)
	Callback for the EventProcessor to take any necessary action when a event input ROOT file is closed.

template<class T >
const T &	getCondition (const std::string &condition_name)
	Access a conditions object for the current event.

TDirectory *	getHistoDirectory ()
	Access/create a directory in the histogram file for this event processor to create histograms and analysis tuples.

void	setStorageHint (framework::StorageControl::Hint hint)
	Mark the current event as having the given storage control hint from this module_.

void	setStorageHint (framework::StorageControl::Hint hint, const std::string &purposeString)
	Mark the current event as having the given storage control hint from this module and the given purpose string.

int	getLogFrequency () const
	Get the current logging frequency from the process.

int	getRunNumber () const
	Get the run number from the process.

std::string	getName () const
	Get the processor name.

void	createHistograms (const std::vector< framework::config::Parameters > &histos)
	Internal function which is used to create histograms passed from the python configuration @parma histos vector of Parameters that configure histograms to create.

Protected Member Functions
ldmx::StraightTrack	seedTracker (const std::tuple< std::array< double, 3 >, ldmx::Measurement, std::optional< ldmx::Measurement > > recoil_one, const std::tuple< std::array< double, 3 >, ldmx::Measurement, std::optional< ldmx::Measurement > > recoil_two, const std::array< double, 3 > ecal_one)

std::tuple< double, double, double, double, std::vector< double > >	fit3DLine (const std::array< double, 3 > &first_recoil, const std::array< double, 3 > &second_recoil, const std::array< double, 3 > &ecal)

double	calculateDistance (const std::array< double, 3 > &point1, const std::array< double, 3 > &point2)

Acts::Vector3	simple3DHitV2 (const ldmx::Measurement &axial, const Acts::Surface &axial_surface, const ldmx::Measurement &stereo, const Acts::Surface &stereo_surface, const ldmx::SimTrackerHit &hitOnTarget, std::vector< ldmx::SimTrackerHit > pair_sim_hits)

std::vector< std::tuple< std::array< double, 3 >, std::tuple< ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit >, std::optional< std::tuple< ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit > > > >	processMeasurements (const std::vector< std::tuple< ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit > > &measurements, const geo::TrackersTrackingGeometry &tg)

double	globalChiSquare (const std::array< double, 3 > &first_sensor, const std::array< double, 3 > &second_sensor, const std::array< double, 3 > &ecal_hit, double a_x, double a_y, double b_x, double b_y)

int	uniqueLayersHit (const std::vector< ldmx::Measurement > &digi_points)

std::array< double, 3 >	convertToLdmxStdArray (const Acts::Vector3 &vec)

std::tuple< Acts::Vector3, Acts::Vector3, Acts::Vector3 >	getSurfaceVectors (const Acts::Surface &surface)

double	dotProduct (const Acts::Vector3 &v1, const Acts::Vector3 &v2)

std::array< double, 3 >	getPointAtZ (std::array< double, 3 > target, std::array< double, 3 > measurement, double z_target)

Protected Member Functions inherited from tracking::reco::TrackingGeometryUser
const Acts::GeometryContext &	geometryContext ()

const Acts::MagneticFieldContext &	magneticFieldContext ()

const Acts::CalibrationContext &	calibrationContext ()

const geo::TrackersTrackingGeometry &	geometry ()

void	loadBField (const std::string &path, const std::vector< double > &map_offset={0., 0., 0.})
	Load the interpolated B-field map from `path` and cache it.

void	loadBField (const std::vector< double > &map_offset={0., 0., 0.})
	Load B-field from the path recorded in the detector GDML.

std::shared_ptr< Acts::MagneticFieldProvider >	bField () const
	Return the loaded B-field provider.

Protected Member Functions inherited from framework::EventProcessor
void	abortEvent ()
	Abort the event immediately.

Protected Attributes
double	processing_time_ {0.}

long	n_events_ {0}

unsigned int	n_seeds_ {0}

std::string	out_seed_collection_ {"LinearRecoilSeedTracks"}
	The name of the output collection of seeds to be stored.

std::string	input_hits_collection_ {"DigiRecoilSimHits"}
	The name of the input hits collection to use in finding seeds..

std::string	input_rec_hits_collection_ {"EcalRecHits"}
	The name of the tagger Tracks (only for Recoil Seeding)

std::string	input_pass_name_ {""}

double	ecal_uncertainty_ {3.87}

double	ecal_distance_threshold_ {10.0}

double	layer12_midpoint_ {12.5}

double	layer23_midpoint_ {20.0}

double	layer34_midpoint_ {27.5}

double	ecal_first_layer_z_threshold_ {250.0}

std::vector< double >	recoil_uncertainty_ {0.006, 0.085}

long	n_missing_ {0}

std::shared_ptr< tracking::sim::TruthMatchingTool >	truth_matching_tool_

Protected Attributes inherited from framework::EventProcessor
HistogramPool	histograms_
	helper object for making and filling histograms

NtupleManager &	ntuple_ {NtupleManager::getInstance()}
	Manager for any ntuples.

logging::logger	the_log_
	The logger for this EventProcessor.

Private Attributes
std::string	next_event_passname_

std::string	sim_particles_passname_

std::string	sim_particles_events_passname_

Detailed Description

Definition at line 26 of file LinearSeedFinder.h.

Constructor & Destructor Documentation

◆ LinearSeedFinder()

tracking::reco::LinearSeedFinder::LinearSeedFinder	(	const std::string &	name,
		framework::Process &	process )

Constructor.

Parameters

name	The name of the instance of this object.
process	The process running this producer.

Definition at line 8 of file LinearSeedFinder.cxx.

10 : TrackingGeometryUser(name, process) {}

framework::Producer::process

virtual void process(Event &event) final

Processing an event for a Producer is calling produce.

Definition EventProcessor.h:278

Member Function Documentation

◆ calculateDistance()

double tracking::reco::LinearSeedFinder::calculateDistance	(	const std::array< double, 3 > &	point1,
		const std::array< double, 3 > &	point2 )

protected

Definition at line 589 of file LinearSeedFinder.cxx.

                                                                          {
  return sqrt(pow(point1[1] - point2[1], 2) + pow(point1[2] - point2[2], 2));
}  // calculateDistance in xy

◆ configure()

void tracking::reco::LinearSeedFinder::configure ( framework::config::Parameters & parameters )

overridevirtual

Configure the processor using the given user specified parameters.

Parameters

parameters Set of parameters used to configure this processor.

Reimplemented from framework::EventProcessor.

Definition at line 16 of file LinearSeedFinder.cxx.

                                                                      {
  // Output seed name
  out_seed_collection_ = parameters.get<std::string>(
      "out_seed_collection", getName() + "LinearRecoilSeedTracks");
 
  // Input strip hits_
  input_hits_collection_ =
      parameters.get<std::string>("input_hits_collection", "DigiRecoilSimHits");
  input_rec_hits_collection_ =
      parameters.get<std::string>("input_rec_hits_collection", "EcalRecHits");
 
  input_pass_name_ = parameters.get<std::string>("input_pass_name", "");
 
  sim_particles_passname_ =
      parameters.get<std::string>("sim_particles_passname");
 
  sim_particles_events_passname_ =
      parameters.get<std::string>("sim_particles_events_passname");
 
  // the uncertainty is sigma_x = 6 microns and sigma_y = 20./sqrt(12)
  recoil_uncertainty_ =
      parameters.get<std::vector<double>>("recoil_uncertainty", {0.006, 0.085});
  ecal_uncertainty_ = parameters.get<double>("ecal_uncertainty", {3.87});
  ecal_distance_threshold_ = parameters.get<double>("ecal_distance_threshold");
  ecal_first_layer_z_threshold_ =
      parameters.get<double>("ecal_first_layer_z_threshold");
 
  layer12_midpoint_ = parameters.get<double>("layer12_midpoint");
  layer23_midpoint_ = parameters.get<double>("layer23_midpoint");
  layer34_midpoint_ = parameters.get<double>("layer34_midpoint");
}

References framework::config::Parameters::get(), framework::EventProcessor::getName(), input_hits_collection_, input_rec_hits_collection_, and out_seed_collection_.

◆ convertToLdmxStdArray()

std::array< double, 3 > tracking::reco::LinearSeedFinder::convertToLdmxStdArray ( const Acts::Vector3 & vec )

protected

Definition at line 428 of file LinearSeedFinder.cxx.

                            {
  return {vec.x(), vec.y(), vec.z()};
}

◆ dotProduct()

double tracking::reco::LinearSeedFinder::dotProduct	(	const Acts::Vector3 &	v1,
		const Acts::Vector3 &	v2 )

protected

Definition at line 536 of file LinearSeedFinder.cxx.

                                                           {
  return v1.dot(v2);
}

◆ fit3DLine()

std::tuple< double, double, double, double, std::vector< double > > tracking::reco::LinearSeedFinder::fit3DLine	(	const std::array< double, 3 > &	first_recoil,
		const std::array< double, 3 > &	second_recoil,
		const std::array< double, 3 > &	ecal )

protected

Definition at line 542 of file LinearSeedFinder.cxx.

                                                             {
  double z_pos1 = first_recoil[0], x_pos1 = first_recoil[1],
         y_pos1 = first_recoil[2];
  double z_pos2 = second_recoil[0], x_pos2 = second_recoil[1],
         y_pos2 = second_recoil[2];
  double z_pos3 = ecal[0], x_pos3 = ecal[1], y_pos3 = ecal[2];
 
  std::array<double, 6> weights = {
      1 / pow(recoil_uncertainty_[0], 2), 1 / pow(recoil_uncertainty_[1], 2),
      1 / pow(recoil_uncertainty_[0], 2), 1 / pow(recoil_uncertainty_[1], 2),
      1 / pow(ecal_uncertainty_, 2),      1 / pow(ecal_uncertainty_, 2)};
 
  Eigen::Matrix<double, 6, 4> a_mat;
  Eigen::Matrix<double, 6, 1> d_vec, w_vec;
 
  // Fill the A matrix (z, 1, 0, 0) for x and (0, 0, z, 1) for y
  a_mat << z_pos1, 1, 0, 0, 0, 0, z_pos1, 1, z_pos2, 1, 0, 0, 0, 0, z_pos2, 1,
      z_pos3, 1, 0, 0, 0, 0, z_pos3, 1;
 
  // Fill the d vector with x and y values
  d_vec << x_pos1, y_pos1, x_pos2, y_pos2, x_pos3, y_pos3;
 
  // Fill the weights vector
  w_vec = Eigen::Matrix<double, 6, 1>(weights.data());
 
  // Solve the weighted least squares system
  Eigen::MatrixXd at_w_a = a_mat.transpose() * w_vec.asDiagonal() * a_mat;
  Eigen::MatrixXd at_w_d = a_mat.transpose() * w_vec.asDiagonal() * d_vec;
  Eigen::VectorXd param_vec = at_w_a.ldlt().solve(at_w_d);
 
  Eigen::Matrix4d covariance_matrix = at_w_a.inverse();
 
  // Store only the upper triangular part of the covariance matrix since it is
  // symmetric
  std::vector<double> covariance_vector = {
      covariance_matrix(0, 0), covariance_matrix(0, 1), covariance_matrix(0, 2),
      covariance_matrix(0, 3), covariance_matrix(1, 1), covariance_matrix(1, 2),
      covariance_matrix(1, 3), covariance_matrix(2, 2), covariance_matrix(2, 3),
      covariance_matrix(3, 3)};
 
  // return {slope_x, intercept_x, slope_y, intercept_y, covariance}
  return {param_vec(0), param_vec(1), param_vec(2), param_vec(3),
          covariance_vector};
}  // fit3DLine

◆ getPointAtZ()

std::array< double, 3 > tracking::reco::LinearSeedFinder::getPointAtZ	(	std::array< double, 3 >	target,
		std::array< double, 3 >	measurement,
		double	z_target )

protected

Definition at line 324 of file LinearSeedFinder.cxx.

                     {
  double slope_x = (measurement[1] - target[0]) / (measurement[0] - target[2]);
  double slope_y = (measurement[2] - target[1]) / (measurement[0] - target[2]);
 
  double intercept_x = target[0] - slope_x * target[2];
  double intercept_y = target[1] - slope_y * target[2];
 
  double x_target = slope_x * z_target + intercept_x;
  double y_target = slope_y * z_target + intercept_y;
 
  return {z_target, x_target, y_target};
}

◆ getSurfaceVectors()

std::tuple< Acts::Vector3, Acts::Vector3, Acts::Vector3 > tracking::reco::LinearSeedFinder::getSurfaceVectors ( const Acts::Surface & surface )

protected

Definition at line 436 of file LinearSeedFinder.cxx.

                                                              {
  Acts::Vector3 dummy{0., 0., 0.};
  Acts::Vector3 u =
      surface.localToGlobal(geometryContext(), Acts::Vector2(1, 0), dummy) -
      surface.center(geometryContext());
  Acts::Vector3 v =
      surface.localToGlobal(geometryContext(), Acts::Vector2(0, 1), dummy) -
      surface.center(geometryContext());
  Acts::Vector3 w = u.cross(v).normalized();
  return {u.normalized(), v.normalized(), w};
}

◆ globalChiSquare()

double tracking::reco::LinearSeedFinder::globalChiSquare	(	const std::array< double, 3 > &	first_sensor,
		const std::array< double, 3 > &	second_sensor,
		const std::array< double, 3 > &	ecal_hit,
		double	a_x,
		double	a_y,
		double	b_x,
		double	b_y )

protected

Definition at line 594 of file LinearSeedFinder.cxx.

                {
  double chi2_x = 0, chi2_y = 0;
  chi2_x += pow(
      (m_x * first_sensor[0] + b_x - first_sensor[1]) / recoil_uncertainty_[0],
      2);
  chi2_y += pow(
      (m_y * first_sensor[0] + b_y - first_sensor[2]) / recoil_uncertainty_[1],
      2);
 
  chi2_x += pow((m_x * second_sensor[0] + b_x - second_sensor[1]) /
                    recoil_uncertainty_[0],
                2);
  chi2_y += pow((m_y * second_sensor[0] + b_y - second_sensor[2]) /
                    recoil_uncertainty_[1],
                2);
 
  chi2_x += pow((m_x * ecal_hit[0] + b_x - ecal_hit[1]) / ecal_uncertainty_, 2);
  chi2_y += pow((m_y * ecal_hit[0] + b_y - ecal_hit[2]) / ecal_uncertainty_, 2);
 
  return chi2_x + chi2_y;
}  // globalChiSquare

◆ onProcessEnd()

void tracking::reco::LinearSeedFinder::onProcessEnd ( )

overridevirtual

Output event statistics.

Reimplemented from framework::EventProcessor.

Definition at line 316 of file LinearSeedFinder.cxx.

                                    {
  ldmx_log(info) << "AVG Time/Event: " << std::fixed << std::setprecision(1)
                 << processing_time_ / n_events_ << " ms";
  ldmx_log(info) << "Total Seeds/Events: " << n_seeds_ << "/" << n_events_;
  ldmx_log(info) << "not enough seed points " << n_missing_;
 
}  // onProcessEnd

◆ onProcessStart()

void tracking::reco::LinearSeedFinder::onProcessStart ( )

overridevirtual

Setup the truth matching.

Reimplemented from framework::EventProcessor.

Definition at line 12 of file LinearSeedFinder.cxx.

                                      {
  truth_matching_tool_ = std::make_shared<tracking::sim::TruthMatchingTool>();
}

◆ processMeasurements()

std::vector< std::tuple< std::array< double, 3 >, std::tuple< ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit >, std::optional< std::tuple< ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit > > > > tracking::reco::LinearSeedFinder::processMeasurements	(	const std::vector< std::tuple< ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit > > &	measurements,
		const geo::TrackersTrackingGeometry &	tg )

protected

Definition at line 344 of file LinearSeedFinder.cxx.

                                           {
  std::vector<
      std::tuple<ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit>>
      axial_measurements;
  std::vector<
      std::tuple<ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit>>
      stereo_measurements;
  std::vector<std::tuple<
      std::array<double, 3>,
      std::tuple<ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit>,
      std::optional<std::tuple<ldmx::Measurement, ldmx::SimTrackerHit,
                               ldmx::SimTrackerHit>>>>
      points_with_measurement;
  Acts::Vector3 dummy{0., 0., 0.};
 
  // Separate measurements into axial and stereo based on layerID
  for (const auto& [measurement, sim_hit, scoring_hit] : measurements) {
    if (measurement.getLayerID() % 2 == 0) {
      axial_measurements.emplace_back(measurement, sim_hit, scoring_hit);
    } else {
      stereo_measurements.emplace_back(measurement, sim_hit, scoring_hit);
    }
  }
 
  // If there are measurements in both axial and stereo layer_:
  // Iterate over axial and stereo measurements to compute 3D points
  if (!axial_measurements.empty() && !stereo_measurements.empty()) {
    for (const auto& axial : axial_measurements) {
      const auto& [axial_meas, axial_hit, axial_sp] = axial;
 
      for (const auto& stereo : stereo_measurements) {
        const auto& [stereo_meas, stereo_hit, stereo_sp] = stereo;
 
        const Acts::Surface* axial_surface =
            tg.getSurface(axial_meas.getLayerID());
        const Acts::Surface* stereo_surface =
            tg.getSurface(stereo_meas.getLayerID());
 
        if (!axial_surface || !stereo_surface) continue;
 
        std::vector<ldmx::SimTrackerHit> sim_hits = {axial_hit, stereo_hit};
        Acts::Vector3 space_point =
            simple3DHitV2(axial_meas, *axial_surface, stereo_meas,
                          *stereo_surface, axial_sp, sim_hits);
 
        points_with_measurement.push_back(
            {convertToLdmxStdArray(space_point), axial, stereo});
      }
    }
  } else if (!axial_measurements.empty()) {
    // if there are only axial measurements, take them to be our sensor points
    for (const auto& axial : axial_measurements) {
      const auto& [axial_meas, axial_hit, axial_sp] = axial;
      const Acts::Surface* axial_surface =
          tg.getSurface(axial_meas.getLayerID());
 
      Acts::Vector3 axial_meas_hit = axial_surface->localToGlobal(
          geometryContext(),
          Acts::Vector2(axial_meas.getLocalPosition()[0], 0.0), dummy);
 
      points_with_measurement.push_back(
          {convertToLdmxStdArray(axial_meas_hit), axial, std::nullopt});
    }
  } else if (!stereo_measurements.empty()) {
    for (const auto& stereo : stereo_measurements) {
      const auto& [stereo_meas, stereo_hit, stereo_sp] = stereo;
      const Acts::Surface* stereo_surface =
          tg.getSurface(stereo_meas.getLayerID());
 
      Acts::Vector3 stereo_meas_hit = stereo_surface->localToGlobal(
          geometryContext(),
          Acts::Vector2(stereo_meas.getLocalPosition()[0], 0.0), dummy);
 
      points_with_measurement.push_back(
          {convertToLdmxStdArray(stereo_meas_hit), stereo, std::nullopt});
    }
  }
  return points_with_measurement;
}

◆ produce()

void tracking::reco::LinearSeedFinder::produce ( framework::Event & event )

overridevirtual

Run the processor and create a collection of results which indicate if a charge particle can be found by the recoil tracker.

Parameters

event The event to process.

Implements framework::Producer.

Definition at line 48 of file LinearSeedFinder.cxx.

                                                    {
  auto start = std::chrono::high_resolution_clock::now();
  std::vector<ldmx::StraightTrack> straight_seed_tracks;
  n_events_++;
  auto tg{geometry()};
 
  const auto& recoil_hits = event.getCollection<ldmx::Measurement>(
      input_hits_collection_, input_pass_name_);
  const auto& ecal_rec_hit = event.getCollection<ldmx::EcalHit>(
      input_rec_hits_collection_, input_pass_name_);
 
  std::vector<std::array<double, 3>> first_layer_ecal_rec_hits;
 
  // Find RecHits at first layer_ of ECal
  for (const auto& x_ecal : ecal_rec_hit) {
    if (x_ecal.getZPos() < ecal_first_layer_z_threshold_) {
      first_layer_ecal_rec_hits.push_back(
          {x_ecal.getZPos(), x_ecal.getXPos(), x_ecal.getYPos()});
    }  // if first layer_ of Ecal
  }  // for positions in ecalRecHit
 
  // Check if we would fit empty seeds, if so: end tracking
  if ((recoil_hits.size() < 2) || (first_layer_ecal_rec_hits.empty()) ||
      (uniqueLayersHit(recoil_hits) < 2)) {
    n_missing_++;
    n_seeds_ += straight_seed_tracks.size();
    event.add(out_seed_collection_, straight_seed_tracks);
    return;
  }
 
  // Setup truth map
  std::map<int, ldmx::SimParticle> particle_map;
  if (event.exists("SimParticles", sim_particles_events_passname_)) {
    particle_map = event.getMap<int, ldmx::SimParticle>(
        "SimParticles", sim_particles_passname_);
    truth_matching_tool_->setup(particle_map, recoil_hits);
  }
 
  std::vector<
      std::tuple<ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit>>
      first_two_layers;
  std::vector<
      std::tuple<ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit>>
      second_two_layers;
 
  const auto& recoil_sim_hits = event.getCollection<ldmx::SimTrackerHit>(
      "RecoilSimHits", input_pass_name_);
  const auto& scoring_hits = event.getCollection<ldmx::SimTrackerHit>(
      "TargetScoringPlaneHits", input_pass_name_);
 
  // Index all sim hits_ by track ID
  std::unordered_map<int, std::vector<const ldmx::SimTrackerHit*>>
      sim_hits_by_track_id;
  for (const auto& hit : recoil_sim_hits) {
    sim_hits_by_track_id[hit.getTrackID()].push_back(&hit);
  }  // for sim hits_
 
  // Index scoring hits_ by track ID (one positive scoring plane hit / ID)
  std::unordered_map<int, const ldmx::SimTrackerHit*> scoring_hit_map;
  for (const auto& sp_hit : scoring_hits) {
    if (sp_hit.getPosition()[2] > 0)
      scoring_hit_map[sp_hit.getTrackID()] = &sp_hit;
  }  // for sp hits_
 
  for (const auto& point : recoil_hits) {
    // x is in tracking coordinates, z is in ldmx coordinates
    float x = point.getGlobalPosition()[0];
    // need to do a size check here since getTrackIds is a std::vector
    // which could be empty
    auto track_ids = point.getTrackIds();
    int track_id = (track_ids.size() > 0) ? track_ids.at(0) : -1;
 
    // get the key value = track_id
    auto sim_range_it = sim_hits_by_track_id.find(track_id);
    if (sim_range_it == sim_hits_by_track_id.end()) continue;
 
    // access map value at track_id
    const auto& sim_hits = sim_range_it->second;
 
    for (const auto* sim_hit : sim_hits) {
      float z = sim_hit->getPosition()[2];
 
      if (x < layer12_midpoint_) {
        if (z < layer12_midpoint_) {
          auto sp_it = scoring_hit_map.find(track_id);
          if (sp_it != scoring_hit_map.end()) {
            first_two_layers.emplace_back(point, *sim_hit, *sp_it->second);
            break;
          }  // add the associated scoring plane hit (will be needed for 3D
             // reconstruction)
        }  // associate 1st layer_ sim hit
      }  // check if recoil hit is 1st layer_
      else if (x < layer23_midpoint_) {
        if (z > layer12_midpoint_ && z < layer23_midpoint_) {
          first_two_layers.emplace_back(point, *sim_hit, ldmx::SimTrackerHit());
          break;
        }  // associate 2nd layer_ sim hit
      }  // check if recoil hit is 2nd layer_
      else if (x < layer34_midpoint_) {
        if (z > layer23_midpoint_ && z < layer34_midpoint_) {
          auto sp_it = scoring_hit_map.find(track_id);
          if (sp_it != scoring_hit_map.end()) {
            second_two_layers.emplace_back(point, *sim_hit, *sp_it->second);
            break;
          }  // add the associated scoring plane hit (will be needed for 3D
             // reconstruction)
        }  // associate 3rd layer_ sim hits_
      }  // check if recoil hit is 3rd layer_
      else {
        if (z > layer34_midpoint_) {
          second_two_layers.emplace_back(point, *sim_hit,
                                         ldmx::SimTrackerHit());
          break;
        }  // associate 4th layer_ sim hits_
      }  // check if recoil hits_ is 4th layer_
 
    }  // loop through simhits
  }  // loop through recoil hits_
 
  // Reconstruct 3D sensor points on which to do fitting
  auto first_sensor_combos = processMeasurements(first_two_layers, tg);
  auto second_sensor_combos = processMeasurements(second_two_layers, tg);
 
  for (const auto& [first_combo_3d_point, first_layer_one, first_layer_two] :
       first_sensor_combos) {
    std::tuple<std::array<double, 3>, ldmx::Measurement,
               std::optional<ldmx::Measurement>>
        first_sensor_point;
 
    if (first_layer_two.has_value()) {
      first_sensor_point = {first_combo_3d_point,
                            std::get<ldmx::Measurement>(first_layer_one),
                            std::get<ldmx::Measurement>(*first_layer_two)};
    }  // check whether we did reconstruction or...
    else {
      first_sensor_point = {first_combo_3d_point,
                            std::get<ldmx::Measurement>(first_layer_one),
                            std::nullopt};
    }  //...we are taking only one layer_ as the measurement (axial or stereo)
 
    for (const auto& [second_combo_3d_point, second_layer_one,
                      second_layer_two] : second_sensor_combos) {
      std::tuple<std::array<double, 3>, ldmx::Measurement,
                 std::optional<ldmx::Measurement>>
          second_sensor_point;
      if (second_layer_two.has_value()) {
        second_sensor_point = {second_combo_3d_point,
                               std::get<ldmx::Measurement>(second_layer_one),
                               std::get<ldmx::Measurement>(*second_layer_two)};
      }  // check whether we did reconstruction or...
      else {
        second_sensor_point = {second_combo_3d_point,
                               std::get<ldmx::Measurement>(second_layer_one),
                               std::nullopt};
      }  //...we are taking only one layer_ as the measurement (axial or stereo)
 
      for (const auto& rec_hit : first_layer_ecal_rec_hits) {
        // Do fitting on 2 sensor + 1 recHit combinations = 1 degree of freedom
        // for linear fit
        ldmx::StraightTrack seed_track =
            seedTracker(first_sensor_point, second_sensor_point, rec_hit);
 
        // Seed passed RecHit distance check, add it
        if (seed_track.getChi2() > 0.0) {
          straight_seed_tracks.push_back(seed_track);
        }  // if chi2 > 0
      }  // for rec_hits
    }  // for second recoil tracker
  }  // for first recoil tracker
 
  n_seeds_ += straight_seed_tracks.size();
  event.add(out_seed_collection_, straight_seed_tracks);
 
  auto end = std::chrono::high_resolution_clock::now();
 
  auto diff = end - start;
  processing_time_ += std::chrono::duration<double, std::milli>(diff).count();
 
  first_layer_ecal_rec_hits.clear();
  straight_seed_tracks.clear();
 
}  // produce

References framework::Event::exists(), input_hits_collection_, input_rec_hits_collection_, and out_seed_collection_.

◆ seedTracker()

ldmx::StraightTrack tracking::reco::LinearSeedFinder::seedTracker	(	const std::tuple< std::array< double, 3 >, ldmx::Measurement, std::optional< ldmx::Measurement > >	recoil_one,
		const std::tuple< std::array< double, 3 >, ldmx::Measurement, std::optional< ldmx::Measurement > >	recoil_two,
		const std::array< double, 3 >	ecal_one )

protected

Definition at line 231 of file LinearSeedFinder.cxx.

                                        {
  auto [sensor1, layer1, layer2] = recoil_one;
  auto [sensor2, layer3, layer4] = recoil_two;
  std::vector<ldmx::Measurement> all_points;
 
  // TODO: in the case where we don't have all 4 hits_, we will be fitting a
  // sensor (weighted average of two layers) + single layer_
  // TODO: or fitting two single layers. Currently, the single layer_ point has
  // the uncertainty of a sensor assigned to it,
  // TODO: but this is not a realistic uncertainty for a single layer_...
  // IF all layers are well-defined, this sequence will add layer1, 2, 3, 4 to
  // the allPoints vector
  all_points.push_back(layer1);
 
  // if layer2 doesn't exist (has_value == False), then the layer1 we added
  // is either layer1 or 2, depending on which one has_value
  if (layer2.has_value()) {
    all_points.push_back(*layer2);
  }
 
  all_points.push_back(layer3);
 
  // if layer4 doesn't exist (has_value == False), then the layer3 we added
  // is either layer3 or 4, depending on which one has_value
  if (layer4.has_value()) {
    all_points.push_back(*layer4);
  }
 
  // Fit the 3 points to a 3D straight line, find track location at first layer_
  // of Ecal, check distance to recHit used in fitting
  // m = slope ; b = intercept
  auto [m_x, b_x, m_y, b_y, seed_cov] = fit3DLine(sensor1, sensor2, ecal_one);
  std::array<double, 3> temp_extrapolated_point = {
      ecal_one[0], m_x * ecal_one[0] + b_x, m_y * ecal_one[0] + b_y};
  double temp_distance = calculateDistance(temp_extrapolated_point, ecal_one);
 
  ldmx::StraightTrack trk = ldmx::StraightTrack();
 
  if (temp_distance < ecal_distance_threshold_) {
    trk.setSlopeX(m_x);
    trk.setInterceptX(b_x);
    trk.setSlopeY(m_y);
    trk.setInterceptY(b_y);
    trk.setTheta(std::atan2(m_y, std::sqrt(1 + m_x * m_x)));
    trk.setPhi(std::atan2(m_x, 1.0));
 
    trk.setAllSensorPoints(all_points);
    trk.setFirstSensorPosition(sensor1);
    trk.setSecondSensorPosition(sensor2);
    trk.setFirstLayerEcalRecHit(ecal_one);
    trk.setDistancetoRecHit(temp_distance);
 
    trk.setTargetLocation(0.0, b_x, b_y);
    trk.setEcalLayer1Location(temp_extrapolated_point);
    trk.setChi2(
        globalChiSquare(sensor1, sensor2, ecal_one, m_x, m_y, b_x, b_y));
    trk.setNhits(3);
    trk.setNdf(1);
 
    trk.setCov(seed_cov);
 
    // truth matching
    if (truth_matching_tool_->configured()) {
      auto truth_info = truth_matching_tool_->truthMatch(all_points);
      trk.setTrackID(truth_info.track_id_);
      trk.setPdgID(truth_info.pdg_id_);
      trk.setTruthProb(truth_info.truth_prob_);
    }
 
    return trk;
 
  }  // if (track is close enough to EcalRecHit)
  else {
    trk.setChi2(-1);
    return trk;
  }  // else (does not pass the threshold)
}  // SeedTracker

◆ simple3DHitV2()

Acts::Vector3 tracking::reco::LinearSeedFinder::simple3DHitV2	(	const ldmx::Measurement &	axial,
		const Acts::Surface &	axial_surface,
		const ldmx::Measurement &	stereo,
		const Acts::Surface &	stereo_surface,
		const ldmx::SimTrackerHit &	hitOnTarget,
		std::vector< ldmx::SimTrackerHit >	pair_sim_hits )

protected

Definition at line 457 of file LinearSeedFinder.cxx.

                                                {
  Acts::Vector3 dummy{0., 0., 0.};
  Acts::Vector3 hit_on_target{target_sp.getPosition()[0],
                              target_sp.getPosition()[1],
                              target_sp.getPosition()[2]};  // x,y,z
 
  Acts::Vector3 axial_true_global{pair_sim_hits[0].getPosition()[0],
                                  pair_sim_hits[0].getPosition()[1],
                                  pair_sim_hits[0].getPosition()[2]};
  // stereo_true_global is unused
  Acts::Vector3 stereo_true_global{pair_sim_hits[1].getPosition()[0],
                                   pair_sim_hits[1].getPosition()[1],
                                   pair_sim_hits[1].getPosition()[2]};
 
  Acts::Vector3 simpart_path = axial_true_global - hit_on_target;
  Acts::Vector3 simpart_unit = simpart_path.normalized();
 
  // Get global positions for strip origins  ....  actually these are in
  // tracking coordinates!
  Acts::Vector3 axial_origin = axial_surface.center(geometryContext());
  Acts::Vector3 stereo_origin = stereo_surface.center(geometryContext());
 
  // the tracking-global vector difference between stereo and axial sensor
  // centers
  Acts::Vector3 delta_sensors = stereo_origin - axial_origin;
 
  // calculate the displacement in tracking global x (need to generalize) by
  // going from tracking x=axial to x=stereo
  double dx_proj = (simpart_unit[0] / simpart_unit[2]) *
                   delta_sensors[0];  // this looks weird because simpart_unit
                                      // is in global-global and delta sensors
                                      // is in tracking-global
 
  // Compute unit vectors for both hits_
  auto [axial_u, axial_v, axial_w] = getSurfaceVectors(axial_surface);
  auto [stereo_u, stereo_v, stereo_w] = getSurfaceVectors(stereo_surface);
  double salpha = dotProduct(axial_v, stereo_u);
  double cosalpha = dotProduct(axial_u, stereo_u);
 
  // Get sensor local measured coordinates
  // Get local position components
  auto [axial_u_value, axial_v_value] = axial.getLocalPosition();
  auto [stereo_u_value, stereo_v_value] = stereo.getLocalPosition();
 
  // Manual correction, since v should always be 0 (insensitive direction)
  axial_v_value = 0.0;
  stereo_v_value = 0.0;
 
  // use the dx_proj as the displacement in u of the axial measurement
  // it's axial_u_value - dx_proj because u is in the -x direction
  //  this calculation is in the axial frame
  double v_intercept_useproj =
      (stereo_u_value - (axial_u_value - dx_proj) * cosalpha) / salpha;
  double u_intercept_useproj = axial_u_value - dx_proj;
 
  // convert to tracking global
  Acts::Vector3 axst_global_useproj = axial_surface.localToGlobal(
      geometryContext(),
      Acts::Vector2(u_intercept_useproj, v_intercept_useproj), dummy);
  Acts::Vector3 dummy_stereo_proj = stereo_surface.localToGlobal(
      geometryContext(),
      Acts::Vector2(u_intercept_useproj, v_intercept_useproj), dummy);
 
  // we want the reconstructed hit to be at the z of the stereo layer
  Acts::Vector3 reconstructed_hit{dummy_stereo_proj[0], axst_global_useproj[1],
                                  axst_global_useproj[2]};
 
  ldmx_log(debug) << "The particle projected axst measured position is "
                     "(compare with stereo sim position): "
                  << reconstructed_hit[0] << ", " << reconstructed_hit[1]
                  << ", " << reconstructed_hit[2] << "\n";
 
  return reconstructed_hit;
}

◆ uniqueLayersHit()

int tracking::reco::LinearSeedFinder::uniqueLayersHit ( const std::vector< ldmx::Measurement > & digi_points )

protected

Definition at line 620 of file LinearSeedFinder.cxx.

                                                   {
  std::vector<ldmx::Measurement> sorted_points = digi_points;
 
  // Sort by z position in the Recoil
  std::sort(sorted_points.begin(), sorted_points.end(),
            [](const ldmx::Measurement& meas1, const ldmx::Measurement& meas2) {
              return meas1.getGlobalPosition()[0] <
                     meas2.getGlobalPosition()[0];
            });
 
  // Remove duplicates to ensure we only keep unique z positions
  auto last = std::unique(
      sorted_points.begin(), sorted_points.end(),
      [](const ldmx::Measurement& meas1, const ldmx::Measurement& meas2) {
        return meas1.getGlobalPosition()[0] == meas2.getGlobalPosition()[0];
      });
 
  // return the number of unique layer_ hits_
  return std::distance(sorted_points.begin(), last);
}  // uniqueLayersHit

Member Data Documentation

◆ ecal_distance_threshold_

double tracking::reco::LinearSeedFinder::ecal_distance_threshold_ {10.0}

protected

Definition at line 136 of file LinearSeedFinder.h.

136{10.0};

◆ ecal_first_layer_z_threshold_

double tracking::reco::LinearSeedFinder::ecal_first_layer_z_threshold_ {250.0}

protected

Definition at line 142 of file LinearSeedFinder.h.

142{250.0};

◆ ecal_uncertainty_

double tracking::reco::LinearSeedFinder::ecal_uncertainty_ {3.87}

protected

Definition at line 134 of file LinearSeedFinder.h.

134{3.87};

◆ input_hits_collection_

std::string tracking::reco::LinearSeedFinder::input_hits_collection_ {"DigiRecoilSimHits"}

protected

The name of the input hits collection to use in finding seeds..

Definition at line 129 of file LinearSeedFinder.h.

129{"DigiRecoilSimHits"};

Referenced by configure(), and produce().

◆ input_pass_name_

std::string tracking::reco::LinearSeedFinder::input_pass_name_ {""}

protected

Definition at line 132 of file LinearSeedFinder.h.

132{""};

◆ input_rec_hits_collection_

std::string tracking::reco::LinearSeedFinder::input_rec_hits_collection_ {"EcalRecHits"}

protected

The name of the tagger Tracks (only for Recoil Seeding)

Definition at line 131 of file LinearSeedFinder.h.

131{"EcalRecHits"};

Referenced by configure(), and produce().

◆ layer12_midpoint_

double tracking::reco::LinearSeedFinder::layer12_midpoint_ {12.5}

protected

Definition at line 139 of file LinearSeedFinder.h.

139{12.5};

◆ layer23_midpoint_

double tracking::reco::LinearSeedFinder::layer23_midpoint_ {20.0}

protected

Definition at line 140 of file LinearSeedFinder.h.

140{20.0};

◆ layer34_midpoint_

double tracking::reco::LinearSeedFinder::layer34_midpoint_ {27.5}

protected

Definition at line 141 of file LinearSeedFinder.h.

141{27.5};

◆ n_events_

long tracking::reco::LinearSeedFinder::n_events_ {0}

protected

Definition at line 123 of file LinearSeedFinder.h.

123{0};

◆ n_missing_

long tracking::reco::LinearSeedFinder::n_missing_ {0}

protected

Definition at line 147 of file LinearSeedFinder.h.

147{0};

◆ n_seeds_

unsigned int tracking::reco::LinearSeedFinder::n_seeds_ {0}

protected

Definition at line 124 of file LinearSeedFinder.h.

124{0};

◆ next_event_passname_

std::string tracking::reco::LinearSeedFinder::next_event_passname_

private

Definition at line 154 of file LinearSeedFinder.h.

◆ out_seed_collection_

std::string tracking::reco::LinearSeedFinder::out_seed_collection_ {"LinearRecoilSeedTracks"}

protected

The name of the output collection of seeds to be stored.

Definition at line 127 of file LinearSeedFinder.h.

127{"LinearRecoilSeedTracks"};

Referenced by configure(), and produce().

◆ processing_time_

double tracking::reco::LinearSeedFinder::processing_time_ {0.}

protected

Definition at line 122 of file LinearSeedFinder.h.

122{0.};

◆ recoil_uncertainty_

std::vector<double> tracking::reco::LinearSeedFinder::recoil_uncertainty_ {0.006, 0.085}

protected

Definition at line 144 of file LinearSeedFinder.h.

144{0.006, 0.085};

◆ sim_particles_events_passname_

std::string tracking::reco::LinearSeedFinder::sim_particles_events_passname_

private

Definition at line 156 of file LinearSeedFinder.h.

◆ sim_particles_passname_

std::string tracking::reco::LinearSeedFinder::sim_particles_passname_

private

Definition at line 155 of file LinearSeedFinder.h.

◆ truth_matching_tool_

std::shared_ptr<tracking::sim::TruthMatchingTool> tracking::reco::LinearSeedFinder::truth_matching_tool_

protected

Initial value:

=

nullptr

Definition at line 150 of file LinearSeedFinder.h.

The documentation for this class was generated from the following files:

Tracking/include/Tracking/Reco/LinearSeedFinder.h
Tracking/src/Tracking/Reco/LinearSeedFinder.cxx

Public Member Functions

Protected Member Functions

Protected Attributes

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ LinearSeedFinder()

Member Function Documentation

◆ calculateDistance()

◆ configure()

◆ convertToLdmxStdArray()

◆ dotProduct()

◆ fit3DLine()

◆ getPointAtZ()

◆ getSurfaceVectors()

◆ globalChiSquare()

◆ onProcessEnd()

◆ onProcessStart()

◆ processMeasurements()

◆ produce()

◆ seedTracker()

◆ simple3DHitV2()

◆ uniqueLayersHit()

Member Data Documentation

◆ ecal_distance_threshold_

◆ ecal_first_layer_z_threshold_

◆ ecal_uncertainty_

◆ input_hits_collection_

◆ input_pass_name_

◆ input_rec_hits_collection_

◆ layer12_midpoint_

◆ layer23_midpoint_

◆ layer34_midpoint_

◆ n_events_

◆ n_missing_

◆ n_seeds_

◆ next_event_passname_

◆ out_seed_collection_

◆ processing_time_

◆ recoil_uncertainty_

◆ sim_particles_events_passname_

◆ sim_particles_passname_

◆ truth_matching_tool_