packages/TrackingDiagnostics/TpcV0CandidateTree.cc

0001 #include "TpcV0CandidateTree.h"
0002
0003 #include <tpctrackreco/TpcTrackHelixFitter.h>
0004 #include <tpctrackreco/TpcTrackKalmanFitter.h>
0005 #include <tpctrackreco/Tpc_PolyCluster.h>
0006 #include <tpctrackreco/Tpc_PolyClusterContainer.h>
0007 #include <tpctrackreco/Tpc_PolyTrack.h>
0008 #include <tpctrackreco/Tpc_PolyTrackContainer.h>
0009 #include <tpctrackreco/Tpc_PolyTrackVertexContainer.h>
0010
0011 #include <trackbase/TrkrDefs.h>
0012
0013 #include <g4main/PHG4EventHeader.h>
0014 #include <g4main/PHG4Hit.h>
0015 #include <g4main/PHG4HitContainer.h>
0016 #include <g4main/PHG4Particle.h>
0017 #include <g4main/PHG4TruthInfoContainer.h>
0018 #include <g4main/PHG4VtxPoint.h>
0019
0020 #include <ffaobjects/EventHeader.h>
0021
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023
0024 #include <phfield/PHField.h>
0025 #include <phfield/PHFieldUtility.h>
0026
0027 #include <phool/PHCompositeNode.h>
0028 #include <phool/getClass.h>
0029
0030 #include <TFile.h>
0031 #include <TTree.h>
0032
0033 #include <Eigen/Dense>
0034
0035 #include <algorithm>
0036 #include <cctype>
0037 #include <chrono>
0038 #include <cmath>
0039 #include <iostream>
0040 #include <limits>
0041 #include <map>
0042 #include <numeric>
0043 #include <tuple>
0044 #include <utility>
0045
0046 namespace
0047 {
0048   constexpr double kPi = 3.14159265358979323846;
0049   constexpr double kPionMass = 0.13957039;
0050   constexpr double kProtonMass = 0.938272088;
0051
0052   template <class T>
0053   constexpr T square(const T &value)
0054   {
0055     return value * value;
0056   }
0057
0058   [[maybe_unused]] int sign_to_charge(const double value)
0059   {
0060     if (!std::isfinite(value) || value == 0.0)
0061     {
0062       return 0;
0063     }
0064     return value > 0.0 ? 1 : -1;
0065   }
0066
0067   double normalize_phi(const double phi)
0068   {
0069     return std::atan2(std::sin(phi), std::cos(phi));
0070   }
0071
0072   double unwrap_to_near(double theta, const double reference)
0073   {
0074     while (theta - reference > kPi)
0075     {
0076       theta -= 2.0 * kPi;
0077     }
0078     while (theta - reference < -kPi)
0079     {
0080       theta += 2.0 * kPi;
0081     }
0082     return theta;
0083   }
0084
0085   bool helix_point_at_beam_radius(const TpcTrackHelix &helix,
0086                                   const double beam_radius,
0087                                   const double theta_hint,
0088                                   double &theta,
0089                                   TpcTrackVec3 &position)
0090   {
0091     if (beam_radius <= 0.0 || helix.radius <= 0.0 ||
0092         !std::isfinite(beam_radius) || !std::isfinite(helix.radius))
0093     {
0094       return false;
0095     }
0096
0097     const double center_radius = std::hypot(helix.cx, helix.cy);
0098     if (center_radius <= 1.0e-12)
0099     {
0100       theta = theta_hint;
0101       position = TpcTrackHelixFitter::point(helix, theta);
0102       return TpcTrackHelixFitter::finite(position);
0103     }
0104
0105     const double rhs = (square(beam_radius) - square(center_radius) - square(helix.radius)) /
0106                        (2.0 * helix.radius * center_radius);
0107     if (rhs < -1.0 - 1.0e-9 || rhs > 1.0 + 1.0e-9)
0108     {
0109       return false;
0110     }
0111
0112     const double clamped_rhs = std::clamp(rhs, -1.0, 1.0);
0113     const double center_phi = std::atan2(helix.cy, helix.cx);
0114     const double delta = std::acos(clamped_rhs);
0115     const double theta_a = unwrap_to_near(center_phi + delta, theta_hint);
0116     const double theta_b = unwrap_to_near(center_phi - delta, theta_hint);
0117     theta = (std::abs(theta_a - theta_hint) <= std::abs(theta_b - theta_hint)) ? theta_a : theta_b;
0118     position = TpcTrackHelixFitter::point(helix, theta);
0119     return TpcTrackHelixFitter::finite(position);
0120   }
0121
0122   bool helix_cluster_residual(const TpcTrackHelix &helix,
0123                               const TpcTrackPoint &point,
0124                               double &previous_theta,
0125                               bool &have_previous_theta,
0126                               TpcTrackVec3 &fit_position,
0127                               double &residual_r,
0128                               double &residual_rphi,
0129                               double &residual_z)
0130   {
0131     const TpcTrackVec3 &cluster = point.position;
0132     const double cluster_r = std::hypot(cluster.x, cluster.y);
0133     double theta_hint = std::atan2(cluster.y - helix.cy,
0134                                    cluster.x - helix.cx);
0135     theta_hint = unwrap_to_near(theta_hint,
0136                                 have_previous_theta ? previous_theta
0137                                                     : helix.theta_first);
0138
0139     double theta = theta_hint;
0140     if (!helix_point_at_beam_radius(helix, cluster_r, theta_hint, theta, fit_position))
0141     {
0142       fit_position = TpcTrackHelixFitter::point(helix, theta_hint);
0143       theta = theta_hint;
0144       if (!TpcTrackHelixFitter::finite(fit_position))
0145       {
0146         return false;
0147       }
0148     }
0149
0150     previous_theta = theta;
0151     have_previous_theta = true;
0152
0153     const double fit_r = std::hypot(fit_position.x, fit_position.y);
0154     const double cluster_phi = std::atan2(cluster.y, cluster.x);
0155     const double fit_phi = std::atan2(fit_position.y, fit_position.x);
0156     residual_r = cluster_r - fit_r;
0157     residual_rphi = cluster_r * normalize_phi(cluster_phi - fit_phi);
0158     residual_z = cluster.z - fit_position.z;
0159     return std::isfinite(residual_r) &&
0160            std::isfinite(residual_rphi) &&
0161            std::isfinite(residual_z);
0162   }
0163 }  // namespace
0164
0165 TpcV0CandidateTree::TpcV0CandidateTree(const std::string &name,
0166                                        const std::string &filename)
0167   : SubsysReco(name)
0168   , m_filename(filename)
0169 {
0170   m_kalman_config.bfield_t = m_bfield_t;
0171   m_kalman_config.point_order = m_point_order;
0172 }
0173
0174 void TpcV0CandidateTree::set_primary_vertex(const double x, const double y, const double z)
0175 {
0176   m_fixed_primary_vertex = {x, y, z};
0177 }
0178
0179 bool TpcV0CandidateTree::set_point_order(const std::string &mode)
0180 {
0181   PointOrder order = PointOrder::Path;
0182   if (!parse_point_order(mode, order))
0183   {
0184     std::cerr << PHWHERE << Name() << ": unknown point order mode '" << mode
0185               << "'. Valid modes are path, input, radius, theta-z, auto." << std::endl;
0186     return false;
0187   }
0188
0189   m_point_order = order;
0190   m_kalman_config.point_order = order;
0191   return true;
0192 }
0193
0194 bool TpcV0CandidateTree::set_track_fit_method(const std::string &mode)
0195 {
0196   std::string lowered = mode;
0197   std::transform(lowered.begin(), lowered.end(), lowered.begin(),
0198                  [](const unsigned char ch)
0199                  { return static_cast<char>(std::tolower(ch)); });
0200
0201   if (lowered == "helix" || lowered == "circle")
0202   {
0203     set_fit_helix(true);
0204     return true;
0205   }
0206   if (lowered == "kalman" || lowered == "kf")
0207   {
0208     set_fit_kalman(true);
0209     return true;
0210   }
0211   if (lowered == "none" || lowered == "line")
0212   {
0213     m_fit_helix_tracks = false;
0214     m_fit_kalman_tracks = false;
0215     return true;
0216   }
0217
0218   std::cerr << PHWHERE << Name() << ": unknown track fit method '" << mode
0219             << "'. Valid methods are helix, kalman, and line." << std::endl;
0220   return false;
0221 }
0222
0223 int TpcV0CandidateTree::Init(PHCompositeNode *topNode)
0224 {
0225   if (m_use_kalman_field_map && m_kalman_config.magnetic_field == nullptr && topNode != nullptr)
0226   {
0227     m_kalman_config.magnetic_field =
0228         findNode::getClass<PHField>(topNode, PHFieldUtility::GetDSTFieldMapNodeName());
0229   }
0230
0231   m_file = new TFile(m_filename.c_str(), "RECREATE");
0232   if (!m_file || m_file->IsZombie())
0233   {
0234     std::cout << Name() << ": failed to create output file " << m_filename << std::endl;
0235     return Fun4AllReturnCodes::ABORTRUN;
0236   }
0237
0238   m_pair_tree = new TTree("pairTree", "TPC truth-point V0 candidates");
0239   m_track_tree = new TTree("trackTree", "TPC track QA before V0 pairing preselection");
0240   if (m_write_cluster_residual_tree)
0241   {
0242     m_cluster_residual_tree = new TTree("clusterResidualTree", "TPC per-cluster residual QA");
0243   }
0244   create_branches();
0245
0246   return Fun4AllReturnCodes::EVENT_OK;
0247 }
0248
0249 int TpcV0CandidateTree::process_event(PHCompositeNode *topNode)
0250 {
0251   const auto event_start = std::chrono::steady_clock::now();
0252   const double kalman_fit_before = m_timing_kalman_fit_seconds;
0253   const double kalman_pca_before = m_timing_kalman_pca_seconds;
0254   const std::uint64_t rkn_propagations_before = m_timing_rkn_propagations;
0255   const std::uint64_t rkn_steps_before = m_timing_rkn_accepted_steps;
0256   const std::uint64_t rkn_retries_before = m_timing_rkn_rejected_trials;
0257   const std::uint64_t rkn_failures_before = m_timing_rkn_failures;
0258   const int run_number = get_run_number(topNode);
0259   const int event_number = get_event_number(topNode);
0260   Vec3 primary_vertex = m_fixed_primary_vertex;
0261   Tpc_PolyTrackVertexContainer *pattern_vertices = nullptr;
0262   std::map<int, Tracklet> tracklet_map;
0263
0264   const auto track_build_start = std::chrono::steady_clock::now();
0265   if (m_use_pattern_cluster_tracks)
0266   {
0267     auto *clusters = findNode::getClass<Tpc_PolyClusterContainer>(
0268         topNode, m_tpc_sa_cluster_node);
0269     auto *tracks = findNode::getClass<Tpc_PolyTrackContainer>(
0270         topNode, m_tpc_sa_track_node);
0271     pattern_vertices = findNode::getClass<Tpc_PolyTrackVertexContainer>(
0272         topNode, m_tpc_sa_track_vertex_node);
0273
0274     if (!clusters || !tracks)
0275     {
0276       if (Verbosity() > 0)
0277       {
0278         std::cout << PHWHERE << Name() << ": missing pattern-reco nodes "
0279                   << m_tpc_sa_cluster_node << "/" << m_tpc_sa_track_node << std::endl;
0280       }
0281       return Fun4AllReturnCodes::EVENT_OK;
0282     }
0283
0284     if (!pattern_vertices && Verbosity() > 1)
0285     {
0286       std::cout << PHWHERE << Name() << ": missing pattern-reco vertex node "
0287                 << m_tpc_sa_track_vertex_node
0288                 << "; trackTree vertex_z will use the configured fallback vertex" << std::endl;
0289     }
0290     tracklet_map = build_pattern_tracklets(clusters, tracks);
0291   }
0292   else
0293   {
0294     auto *truth_points = findNode::getClass<PHG4HitContainer>(topNode, m_truth_point_node);
0295     if (!truth_points)
0296     {
0297       if (Verbosity() > 0)
0298       {
0299         std::cout << PHWHERE << Name() << ": missing truth point node "
0300                   << m_truth_point_node << std::endl;
0301       }
0302       return Fun4AllReturnCodes::EVENT_OK;
0303     }
0304
0305     auto *truth_info = findNode::getClass<PHG4TruthInfoContainer>(topNode, m_truth_info_node);
0306     if (!truth_info && Verbosity() > 0)
0307     {
0308       std::cout << PHWHERE << Name() << ": missing truth info node "
0309                 << m_truth_info_node << "; no charged truth tracklets can be built" << std::endl;
0310     }
0311
0312     primary_vertex = get_primary_vertex(truth_info);
0313     tracklet_map = build_tracklets(truth_points, truth_info);
0314   }
0315   const double track_build_seconds = std::chrono::duration<double>(
0316                                          std::chrono::steady_clock::now() - track_build_start)
0317                                          .count();
0318   if (m_print_timing)
0319   {
0320     std::cout << "[V0TimingStage] run=" << run_number
0321               << " event=" << event_number
0322               << " stage=track_build_done"
0323               << " tracks=" << tracklet_map.size()
0324               << " stage_s=" << track_build_seconds
0325               << " kalman_fit_s=" << (m_timing_kalman_fit_seconds - kalman_fit_before)
0326               << " rkn_propagations=" << (m_timing_rkn_propagations - rkn_propagations_before)
0327               << " rkn_steps=" << (m_timing_rkn_accepted_steps - rkn_steps_before)
0328               << std::endl;
0329   }
0330
0331   const auto dca_cache_start = std::chrono::steady_clock::now();
0332   for (auto &entry : tracklet_map)
0333   {
0334     auto &tracklet = entry.second;
0335     tracklet.has_beamline_pca = track_pca_to_xy(
0336         tracklet, Vec3{0.0, 0.0, 0.0}, tracklet.beamline_pca, tracklet.rdca_zero);
0337     tracklet.has_pattern_vertex = choose_pattern_collision_vertex(
0338         tracklet, pattern_vertices, tracklet.pattern_vertex,
0339         tracklet.pattern_vertex_z_rms, tracklet.pattern_vertex_ntracks);
0340     const Vec3 &dca_vertex = tracklet.has_pattern_vertex
0341                                  ? tracklet.pattern_vertex
0342                                  : primary_vertex;
0343     tracklet.vertex_dca = fitted_track_dca_to_vertex(tracklet, dca_vertex);
0344     tracklet.has_vertex_dca =
0345         std::isfinite(tracklet.vertex_dca.first) &&
0346         std::isfinite(tracklet.vertex_dca.second);
0347   }
0348   const double dca_cache_seconds = std::chrono::duration<double>(
0349                                        std::chrono::steady_clock::now() - dca_cache_start)
0350                                        .count();
0351   m_timing_dca_cache_seconds += dca_cache_seconds;
0352   if (m_print_timing)
0353   {
0354     std::cout << "[V0TimingStage] run=" << run_number
0355               << " event=" << event_number
0356               << " stage=dca_cache_done"
0357               << " tracks=" << tracklet_map.size()
0358               << " stage_s=" << dca_cache_seconds
0359               << std::endl;
0360   }
0361
0362   std::vector<const Tracklet *> tracklets;
0363   tracklets.reserve(tracklet_map.size());
0364   for (const auto &entry : tracklet_map)
0365   {
0366     tracklets.push_back(&entry.second);
0367   }
0368
0369   const auto track_qa_start = std::chrono::steady_clock::now();
0370   for (const auto *tracklet : tracklets)
0371   {
0372     fill_track_row(*tracklet, primary_vertex, run_number, event_number);
0373     if (m_cluster_residual_tree)
0374     {
0375       fill_cluster_residual_rows(*tracklet, primary_vertex, run_number, event_number);
0376     }
0377   }
0378   const double track_qa_seconds = std::chrono::duration<double>(
0379                                       std::chrono::steady_clock::now() - track_qa_start)
0380                                       .count();
0381   if (m_print_timing)
0382   {
0383     std::cout << "[V0TimingStage] run=" << run_number
0384               << " event=" << event_number
0385               << " stage=track_qa_done"
0386               << " stage_s=" << track_qa_seconds
0387               << std::endl;
0388   }
0389
0390   const auto pair_loop_start = std::chrono::steady_clock::now();
0391   std::uint64_t event_pairs_processed = 0;
0392   const std::uint64_t event_pairs_total =
0393       tracklets.size() > 1
0394           ? static_cast<std::uint64_t>(tracklets.size()) *
0395                 static_cast<std::uint64_t>(tracklets.size() - 1) / 2
0396           : 0;
0397   if (m_print_timing)
0398   {
0399     std::cout << "[V0TimingStage] run=" << run_number
0400               << " event=" << event_number
0401               << " stage=pair_loop_start"
0402               << " pairs=" << event_pairs_total
0403               << std::endl;
0404   }
0405   for (std::size_t i = 0; i < tracklets.size(); ++i)
0406   {
0407     for (std::size_t j = i + 1; j < tracklets.size(); ++j)
0408     {
0409       make_pair_row(*tracklets[i], *tracklets[j], primary_vertex, run_number, event_number);
0410       ++event_pairs_processed;
0411       if (m_print_timing &&
0412           (event_pairs_processed == 1 || event_pairs_processed % 1000 == 0))
0413       {
0414         std::cout << "[V0TimingPair] run=" << run_number
0415                   << " event=" << event_number
0416                   << " done=" << event_pairs_processed
0417                   << " total=" << event_pairs_total
0418                   << " pair_loop_s=" << std::chrono::duration<double>(std::chrono::steady_clock::now() - pair_loop_start).count()
0419                   << " kalman_pca_s=" << (m_timing_kalman_pca_seconds - kalman_pca_before)
0420                   << std::endl;
0421       }
0422     }
0423   }
0424   const double pair_loop_seconds = std::chrono::duration<double>(
0425                                        std::chrono::steady_clock::now() - pair_loop_start)
0426                                        .count();
0427   const double total_seconds = std::chrono::duration<double>(
0428                                    std::chrono::steady_clock::now() - event_start)
0429                                    .count();
0430   const double kalman_fit_seconds = m_timing_kalman_fit_seconds - kalman_fit_before;
0431   const double kalman_pca_seconds = m_timing_kalman_pca_seconds - kalman_pca_before;
0432
0433   ++m_timing_events;
0434   m_timing_total_seconds += total_seconds;
0435   m_timing_track_build_seconds += track_build_seconds;
0436   m_timing_track_qa_seconds += track_qa_seconds;
0437   m_timing_pair_loop_seconds += pair_loop_seconds;
0438
0439   if (m_print_timing)
0440   {
0441     std::cout << "[V0Timing] run=" << run_number
0442               << " event=" << event_number
0443               << " tracks=" << tracklets.size()
0444               << " total_s=" << total_seconds
0445               << " track_build_s=" << track_build_seconds
0446               << " kalman_fit_s=" << kalman_fit_seconds
0447               << " track_qa_s=" << track_qa_seconds
0448               << " pair_loop_s=" << pair_loop_seconds
0449               << " kalman_pca_s=" << kalman_pca_seconds
0450               << " rkn_propagations=" << (m_timing_rkn_propagations - rkn_propagations_before)
0451               << " rkn_steps=" << (m_timing_rkn_accepted_steps - rkn_steps_before)
0452               << " rkn_retries=" << (m_timing_rkn_rejected_trials - rkn_retries_before)
0453               << " rkn_failures=" << (m_timing_rkn_failures - rkn_failures_before)
0454               << std::endl;
0455   }
0456
0457   return Fun4AllReturnCodes::EVENT_OK;
0458 }
0459
0460 int TpcV0CandidateTree::End(PHCompositeNode * /*topNode*/)
0461 {
0462   if (m_file)
0463   {
0464     m_file->cd();
0465     if (m_pair_tree)
0466     {
0467       m_pair_tree->Write();
0468     }
0469     if (m_track_tree)
0470     {
0471       m_track_tree->Write();
0472     }
0473     if (m_cluster_residual_tree)
0474     {
0475       m_cluster_residual_tree->Write();
0476     }
0477     m_file->Close();
0478     delete m_file;
0479     m_file = nullptr;
0480   }
0481
0482   if (Verbosity() > 0)
0483   {
0484     std::cout << Name() << ": pair counters: raw=" << m_counter_raw_pairs
0485               << " reject_charge=" << m_counter_reject_charge
0486               << " reject_preselection=" << m_counter_reject_preselection
0487               << " reject_pca=" << m_counter_reject_pca
0488               << " reject_pointing=" << m_counter_reject_pointing
0489               << " reject_ap=" << m_counter_reject_ap
0490               << " reject_pair_selection=" << m_counter_reject_pair_selection
0491               << " written=" << m_counter_written
0492               << " tracks_written=" << m_counter_tracks_written
0493               << " reject_helix_anchor=" << m_counter_reject_helix_anchor
0494               << " cluster_residuals_written=" << m_counter_cluster_residuals_written
0495               << std::endl;
0496   }
0497
0498   if (m_print_timing)
0499   {
0500     std::cout << "[V0TimingSummary] events=" << m_timing_events
0501               << " total_s=" << m_timing_total_seconds
0502               << " track_build_s=" << m_timing_track_build_seconds
0503               << " kalman_fits=" << m_timing_kalman_fits
0504               << " kalman_fit_s=" << m_timing_kalman_fit_seconds
0505               << " rkn_s=" << m_timing_rkn_seconds
0506               << " dca_cache_s=" << m_timing_dca_cache_seconds
0507               << " track_qa_s=" << m_timing_track_qa_seconds
0508               << " pair_loop_s=" << m_timing_pair_loop_seconds
0509               << " kalman_pca_s=" << m_timing_kalman_pca_seconds
0510               << " rkn_propagations=" << m_timing_rkn_propagations
0511               << " rkn_steps=" << m_timing_rkn_accepted_steps
0512               << " rkn_retries=" << m_timing_rkn_rejected_trials
0513               << " rkn_failures=" << m_timing_rkn_failures
0514               << std::endl;
0515   }
0516
0517   return Fun4AllReturnCodes::EVENT_OK;
0518 }
0519
0520 int TpcV0CandidateTree::get_event_number(PHCompositeNode *topNode) const
0521 {
0522   if (auto *event_header = findNode::getClass<EventHeader>(topNode, "EventHeader"))
0523   {
0524     return event_header->get_EvtSequence();
0525   }
0526   if (auto *g4_event_header = findNode::getClass<PHG4EventHeader>(topNode, "EventHeader"))
0527   {
0528     return g4_event_header->get_EvtSequence();
0529   }
0530   return 0;
0531 }
0532
0533 int TpcV0CandidateTree::get_run_number(PHCompositeNode *topNode) const
0534 {
0535   if (auto *event_header = findNode::getClass<EventHeader>(topNode, "EventHeader"))
0536   {
0537     return event_header->get_RunNumber();
0538   }
0539   return 1;
0540 }
0541
0542 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::get_primary_vertex(PHG4TruthInfoContainer *truth_info) const
0543 {
0544   if (!m_use_truth_primary_vertex || !truth_info)
0545   {
0546     return m_fixed_primary_vertex;
0547   }
0548
0549   const auto vtx_range = truth_info->GetPrimaryVtxRange();
0550   for (auto iter = vtx_range.first; iter != vtx_range.second; ++iter)
0551   {
0552     const auto *vtx = iter->second;
0553     if (vtx)
0554     {
0555       return {vtx->get_x(), vtx->get_y(), vtx->get_z()};
0556     }
0557   }
0558
0559   return m_fixed_primary_vertex;
0560 }
0561
0562 std::map<int, TpcV0CandidateTree::Tracklet> TpcV0CandidateTree::build_tracklets(
0563     PHG4HitContainer *truth_points,
0564     PHG4TruthInfoContainer *truth_info) const
0565 {
0566   std::map<int, Tracklet> tracklets;
0567   if (!truth_points)
0568   {
0569     return tracklets;
0570   }
0571
0572   const auto hit_range = truth_points->getHits();
0573   for (auto hit_iter = hit_range.first; hit_iter != hit_range.second; ++hit_iter)
0574   {
0575     const PHG4Hit *hit = hit_iter->second;
0576     if (!hit)
0577     {
0578       continue;
0579     }
0580
0581     const int track_id = hit->get_trkid();
0582     Tracklet &tracklet = tracklets[track_id];
0583     if (tracklet.points.empty())
0584     {
0585       tracklet.track_id = track_id;
0586       tracklet.shower_id = hit->get_shower_id();
0587     }
0588
0589     TruthPoint point;
0590     point.track_id = track_id;
0591     point.shower_id = hit->get_shower_id();
0592     point.layer = static_cast<int>(hit->get_layer());
0593     point.position = {hit->get_x(0), hit->get_y(0), hit->get_z(0)};
0594     point.momentum = {hit->get_px(0), hit->get_py(0), hit->get_pz(0)};
0595     point.t = hit->get_t(0);
0596     point.path = hit->get_path_length();
0597     if (finite(point.position) && finite(point.momentum))
0598     {
0599       tracklet.points.push_back(point);
0600     }
0601   }
0602
0603   for (auto iter = tracklets.begin(); iter != tracklets.end();)
0604   {
0605     Tracklet &tracklet = iter->second;
0606     order_track_points(tracklet.points, m_point_order);
0607     tracklet.npoints = static_cast<int>(tracklet.points.size());
0608     tracklet.ntpc_clusters = static_cast<unsigned int>(tracklet.points.size());
0609     if (!tracklet.points.empty())
0610     {
0611       tracklet.side = tracklet.points.front().position.z < 0.0 ? 0 : 1;
0612     }
0613
0614     if (tracklet.npoints < m_min_points)
0615     {
0616       iter = tracklets.erase(iter);
0617       continue;
0618     }
0619
0620     if (truth_info)
0621     {
0622       const PHG4Particle *particle = truth_info->GetParticle(tracklet.track_id);
0623       if (particle)
0624       {
0625         tracklet.pid = particle->get_pid();
0626         tracklet.parent_id = particle->get_parent_id();
0627         tracklet.primary_id = particle->get_primary_id();
0628         tracklet.vtx_id = particle->get_vtx_id();
0629         tracklet.barcode = particle->get_barcode();
0630         tracklet.embed_id = truth_info->isEmbeded(tracklet.track_id);
0631         tracklet.is_primary = truth_info->is_primary(particle) ? 1 : 0;
0632         tracklet.charge = pdg_charge(tracklet.pid);
0633         tracklet.truth_momentum = {particle->get_px(), particle->get_py(), particle->get_pz()};
0634         tracklet.truth_e = particle->get_e();
0635
0636         if (const PHG4Particle *parent = truth_info->GetParticle(tracklet.parent_id))
0637         {
0638           tracklet.parent_pid = parent->get_pid();
0639         }
0640
0641         if (auto *vtx = truth_info->GetVtx(tracklet.vtx_id))
0642         {
0643           tracklet.truth_vertex = {vtx->get_x(), vtx->get_y(), vtx->get_z()};
0644           tracklet.truth_vt = vtx->get_t();
0645         }
0646       }
0647     }
0648
0649     if (tracklet.charge == 0)
0650     {
0651       iter = tracklets.erase(iter);
0652       continue;
0653     }
0654
0655     if (m_fit_kalman_tracks)
0656     {
0657       tracklet.has_kalman = fit_kalman(tracklet.points, tracklet.charge, tracklet.kalman);
0658       if (!tracklet.has_kalman || tracklet.kalman.states_smoothed.empty())
0659       {
0660         iter = tracklets.erase(iter);
0661         continue;
0662       }
0663
0664       const auto state = TpcTrackKalmanFitter::propagation_state(tracklet.kalman, Vec3{});
0665       tracklet.position = TpcTrackKalmanFitter::state_position(state);
0666       tracklet.momentum = TpcTrackKalmanFitter::state_momentum(state);
0667     }
0668     else if (m_fit_helix_tracks)
0669     {
0670       tracklet.has_helix = fit_helix(tracklet.points, m_fit_first_points,
0671                                      tracklet.charge, m_bfield_t, tracklet.helix);
0672       if (!tracklet.has_helix)
0673       {
0674         iter = tracklets.erase(iter);
0675         continue;
0676       }
0677       tracklet.position = helix_point(tracklet.helix, tracklet.helix.theta_first);
0678       tracklet.momentum = helix_momentum(tracklet.helix, tracklet.helix.theta_first);
0679     }
0680     else
0681     {
0682       tracklet.position = tracklet.points.front().position;
0683       tracklet.momentum = tracklet.points.front().momentum;
0684     }
0685
0686     assign_fit_quality(tracklet);
0687
0688     if (!finite(tracklet.truth_momentum) || norm(tracklet.truth_momentum) <= 0.0)
0689     {
0690       tracklet.truth_momentum = tracklet.momentum;
0691     }
0692
0693     ++iter;
0694   }
0695
0696   return tracklets;
0697 }
0698
0699 std::map<int, TpcV0CandidateTree::Tracklet> TpcV0CandidateTree::build_pattern_tracklets(
0700     Tpc_PolyClusterContainer *clusters,
0701     Tpc_PolyTrackContainer *tracks) const
0702 {
0703   std::map<int, Tracklet> tracklets;
0704
0705   if (!clusters || !tracks)
0706   {
0707     return tracklets;
0708   }
0709
0710   std::map<unsigned int, std::vector<const Tpc_PolyCluster *>> clusters_by_source_id;
0711   for (unsigned int icluster = 0; icluster < clusters->size(); ++icluster)
0712   {
0713     const Tpc_PolyCluster *cluster = clusters->get_cluster(icluster);
0714     if (!cluster || !cluster->isValid())
0715     {
0716       continue;
0717     }
0718     clusters_by_source_id[cluster->get_source_assembled_track_id()].push_back(cluster);
0719   }
0720
0721   for (unsigned int itrack = 0; itrack < tracks->size(); ++itrack)
0722   {
0723     const Tpc_PolyTrack *track = tracks->get_track(itrack);
0724     if (!track || !track->isValid() || track->get_fit_status() == 0)
0725     {
0726       continue;
0727     }
0728
0729     const unsigned int source_id = track->get_source_assembled_track_id();
0730     const auto cluster_iter = clusters_by_source_id.find(source_id);
0731     if (cluster_iter == clusters_by_source_id.end())
0732     {
0733       continue;
0734     }
0735
0736     const unsigned int pattern_track_id = track->get_track_id();
0737     const int track_id = static_cast<int>(pattern_track_id != 0 ? pattern_track_id : itrack + 1);
0738
0739     Tracklet tracklet;
0740     tracklet.track_id = track_id;
0741     tracklet.shower_id = static_cast<int>(source_id);
0742     tracklet.charge = sign_to_charge(track->get_charge());
0743     tracklet.position = {track->get_x(), track->get_y(), track->get_z()};
0744     tracklet.momentum = {track->get_px(), track->get_py(), track->get_pz()};
0745     tracklet.truth_momentum = tracklet.momentum;
0746     tracklet.dedx = track->get_dedx();
0747     tracklet.has_dedx = std::isfinite(tracklet.dedx);
0748
0749     std::map<int, unsigned int> side_counts;
0750     const auto &track_clusters = cluster_iter->second;
0751     tracklet.points.reserve(track_clusters.size());
0752     for (unsigned int icluster = 0; icluster < track_clusters.size(); ++icluster)
0753     {
0754       const Tpc_PolyCluster *cluster = track_clusters[icluster];
0755       if (!cluster || !cluster->isValid())
0756       {
0757         continue;
0758       }
0759
0760       TruthPoint point;
0761       point.track_id = track_id;
0762       point.shower_id = static_cast<int>(source_id);
0763       point.layer = -1;
0764       if (cluster->size_hits() > 0)
0765       {
0766         point.layer = static_cast<int>(TrkrDefs::getLayer(cluster->get_hit_index(0).first));
0767       }
0768       point.position = {cluster->get_centroid_x(),
0769                         cluster->get_centroid_y(),
0770                         cluster->get_centroid_z()};
0771       point.momentum = tracklet.momentum;
0772       point.t = 0.0;
0773       point.path = point.layer >= 0 ? static_cast<double>(point.layer)
0774                                     : static_cast<double>(icluster);
0775       if (finite(point.position))
0776       {
0777         tracklet.points.push_back(point);
0778         ++side_counts[cluster->get_side()];
0779       }
0780     }
0781
0782     if (!side_counts.empty())
0783     {
0784       tracklet.side = std::max_element(
0785                           side_counts.begin(), side_counts.end(),
0786                           [](const auto &lhs, const auto &rhs)
0787                           { return lhs.second < rhs.second; })
0788                           ->first;
0789     }
0790     tracklet.ntpc_clusters = track->get_nclusters() > 0
0791                                  ? track->get_nclusters()
0792                                  : static_cast<unsigned int>(tracklet.points.size());
0793
0794     const bool has_upstream_state = finite(tracklet.position) &&
0795                                     finite(tracklet.momentum) &&
0796                                     norm(tracklet.momentum) > 0.0;
0797     if (!finalize_pattern_tracklet(tracklet, has_upstream_state))
0798     {
0799       continue;
0800     }
0801
0802     tracklet.truth_momentum = tracklet.momentum;
0803     tracklets[track_id] = std::move(tracklet);
0804   }
0805   return tracklets;
0806 }
0807
0808 bool TpcV0CandidateTree::finalize_pattern_tracklet(Tracklet &tracklet,
0809                                                    const bool has_upstream_state) const
0810 {
0811   order_track_points(tracklet.points, m_point_order);
0812   tracklet.npoints = static_cast<int>(tracklet.points.size());
0813   if (tracklet.npoints < m_min_points || tracklet.charge == 0)
0814   {
0815     return false;
0816   }
0817
0818   if (m_fit_kalman_tracks)
0819   {
0820     tracklet.has_kalman = fit_kalman(tracklet.points, tracklet.charge, tracklet.kalman);
0821     if (!tracklet.has_kalman || tracklet.kalman.states_smoothed.empty())
0822     {
0823       return false;
0824     }
0825
0826     const auto state = TpcTrackKalmanFitter::propagation_state(tracklet.kalman, Vec3{});
0827     tracklet.position = TpcTrackKalmanFitter::state_position(state);
0828     tracklet.momentum = TpcTrackKalmanFitter::state_momentum(state);
0829   }
0830   else if (m_use_final_track_helix && has_upstream_state)
0831   {
0832     tracklet.has_helix = helix_from_state(tracklet.position, tracklet.momentum,
0833                                           tracklet.charge, m_bfield_t, tracklet.helix);
0834     if (!tracklet.has_helix)
0835     {
0836       return false;
0837     }
0838
0839     tracklet.has_helix_search_range =
0840         TpcTrackHelixFitter::measurement_anchored_search_range(
0841             tracklet.helix, tracklet.points,
0842             m_final_track_helix_max_upstream_cm,
0843             m_final_track_helix_downstream_margin_cm,
0844             tracklet.helix_search_range);
0845     if (!tracklet.has_helix_search_range)
0846     {
0847       ++m_counter_reject_helix_anchor;
0848       return false;
0849     }
0850
0851     tracklet.position = helix_point(tracklet.helix, tracklet.helix.theta_first);
0852     tracklet.momentum = helix_momentum(tracklet.helix, tracklet.helix.theta_first);
0853   }
0854   else if (m_fit_helix_tracks)
0855   {
0856     tracklet.has_helix = fit_helix(tracklet.points, m_fit_first_points,
0857                                    tracklet.charge, m_bfield_t, tracklet.helix);
0858     if (!tracklet.has_helix)
0859     {
0860       return false;
0861     }
0862
0863     tracklet.position = helix_point(tracklet.helix, tracklet.helix.theta_first);
0864     tracklet.momentum = helix_momentum(tracklet.helix, tracklet.helix.theta_first);
0865   }
0866   else if (!has_upstream_state)
0867   {
0868     if (tracklet.points.size() < 2)
0869     {
0870       return false;
0871     }
0872     tracklet.position = tracklet.points.front().position;
0873     tracklet.momentum = subtract(tracklet.points[1].position,
0874                                  tracklet.points.front().position);
0875   }
0876
0877   assign_fit_quality(tracklet);
0878   return finite(tracklet.position) && finite(tracklet.momentum) &&
0879          norm(tracklet.momentum) > 0.0;
0880 }
0881
0882 bool TpcV0CandidateTree::track_pca_to_xy(const Tracklet &tracklet,
0883                                          const Vec3 &beamline,
0884                                          Vec3 &pca,
0885                                          double &signed_dca_xy) const
0886 {
0887   const double nan = std::numeric_limits<double>::quiet_NaN();
0888   pca = {nan, nan, nan};
0889   signed_dca_xy = nan;
0890
0891   if (tracklet.has_kalman && !tracklet.kalman.states_smoothed.empty())
0892   {
0893     const auto state = TpcTrackKalmanFitter::propagation_state(tracklet.kalman, beamline);
0894     const Vec3 position = TpcTrackKalmanFitter::state_position(state);
0895     const Vec3 momentum = TpcTrackKalmanFitter::state_momentum(state);
0896     const double qop_t = state[TpcTrackKalmanFitter::QOverPt];
0897     const double omega = 0.003 * tracklet.kalman.bfield_t * qop_t;
0898     if (std::abs(omega) >= 1.0e-10)
0899     {
0900       const double radius = std::abs(1.0 / omega);
0901       const double center_x = state[TpcTrackKalmanFitter::X] -
0902                               std::sin(state[TpcTrackKalmanFitter::Phi]) / omega;
0903       const double center_y = state[TpcTrackKalmanFitter::Y] +
0904                               std::cos(state[TpcTrackKalmanFitter::Phi]) / omega;
0905       const double dx = beamline.x - center_x;
0906       const double dy = beamline.y - center_y;
0907       const double center_distance = std::hypot(dx, dy);
0908       if (center_distance > 1.0e-12)
0909       {
0910         const double closest_x = center_x + radius * dx / center_distance;
0911         const double closest_y = center_y + radius * dy / center_distance;
0912         const double theta0 = std::atan2(position.y - center_y, position.x - center_x);
0913         const double theta_closest = std::atan2(closest_y - center_y,
0914                                                 closest_x - center_x);
0915         const double path_cm = normalize_phi(theta_closest - theta0) / omega;
0916         TpcKalmanConfig config = m_kalman_config;
0917         config.bfield_t = tracklet.kalman.bfield_t;
0918         config.magnetic_field = tracklet.kalman.magnetic_field;
0919         config.analytic_uniform_propagation = tracklet.kalman.analytic_uniform_propagation;
0920         const auto closest_state = TpcTrackKalmanFitter::propagate_state(
0921             state, path_cm, config, tracklet.kalman.mass_gev);
0922         pca = TpcTrackKalmanFitter::state_position(closest_state);
0923         signed_dca_xy = center_distance - radius;
0924         return finite(pca) && std::isfinite(signed_dca_xy);
0925       }
0926
0927       pca = position;
0928       signed_dca_xy = -radius;
0929       return finite(pca);
0930     }
0931
0932     const double transverse_momentum2 = square(momentum.x) + square(momentum.y);
0933     if (transverse_momentum2 <= 0.0)
0934     {
0935       return false;
0936     }
0937     const Vec3 relative = subtract(position, beamline);
0938     const double scale_to_pca = -(relative.x * momentum.x + relative.y * momentum.y) /
0939                                 transverse_momentum2;
0940     pca = add(position, scale(momentum, scale_to_pca));
0941     signed_dca_xy = (relative.x * momentum.y - relative.y * momentum.x) /
0942                     std::sqrt(transverse_momentum2);
0943     return finite(pca) && std::isfinite(signed_dca_xy);
0944   }
0945
0946   if (tracklet.has_helix && tracklet.helix.radius > 0.0)
0947   {
0948     const double dx = beamline.x - tracklet.helix.cx;
0949     const double dy = beamline.y - tracklet.helix.cy;
0950     const double center_distance = std::hypot(dx, dy);
0951     double theta = tracklet.helix.theta_first;
0952     if (center_distance > 1.0e-12)
0953     {
0954       const double theta_raw = std::atan2(dy, dx);
0955       theta = unwrap_to_near(theta_raw, tracklet.helix.theta_first);
0956     }
0957     pca = helix_point(tracklet.helix, theta);
0958     signed_dca_xy = center_distance - tracklet.helix.radius;
0959     return finite(pca) && std::isfinite(signed_dca_xy);
0960   }
0961
0962   const double transverse_momentum2 = square(tracklet.momentum.x) + square(tracklet.momentum.y);
0963   if (transverse_momentum2 <= 0.0)
0964   {
0965     return false;
0966   }
0967   const Vec3 relative = subtract(tracklet.position, beamline);
0968   const double scale_to_pca = -(relative.x * tracklet.momentum.x +
0969                                 relative.y * tracklet.momentum.y) /
0970                               transverse_momentum2;
0971   pca = add(tracklet.position, scale(tracklet.momentum, scale_to_pca));
0972   signed_dca_xy = (relative.x * tracklet.momentum.y -
0973                    relative.y * tracklet.momentum.x) /
0974                   std::sqrt(transverse_momentum2);
0975   return finite(pca) && std::isfinite(signed_dca_xy);
0976 }
0977
0978 bool TpcV0CandidateTree::choose_pattern_collision_vertex(
0979     const Tracklet &tracklet,
0980     Tpc_PolyTrackVertexContainer *vertices,
0981     Vec3 &vertex,
0982     double &z_rms,
0983     unsigned int &ntracks) const
0984 {
0985   if (!vertices || vertices->get_collision_vertex_valid() == 0)
0986   {
0987     return false;
0988   }
0989
0990   double best_dz = std::numeric_limits<double>::max();
0991   const unsigned int count = vertices->get_collision_vertex_count();
0992   for (unsigned int index = 0; index < count; ++index)
0993   {
0994     const Vec3 candidate{vertices->get_collision_x(index),
0995                          vertices->get_collision_y(index),
0996                          vertices->get_collision_z(index)};
0997     if (!finite(candidate))
0998     {
0999       continue;
1000     }
1001
1002     Vec3 candidate_pca;
1003     double candidate_dca = 0.0;
1004     if (!track_pca_to_xy(tracklet, candidate, candidate_pca, candidate_dca))
1005     {
1006       continue;
1007     }
1008     const double dz = std::abs(candidate_pca.z - candidate.z);
1009     if (dz < best_dz)
1010     {
1011       best_dz = dz;
1012       vertex = candidate;
1013       z_rms = vertices->get_collision_z_rms(index);
1014       ntracks = vertices->get_collision_ntracks(index);
1015     }
1016   }
1017   return best_dz < std::numeric_limits<double>::max();
1018 }
1019
1020 bool TpcV0CandidateTree::make_pair_row(const Tracklet &track1, const Tracklet &track2,
1021                                        const Vec3 &primary_vertex,
1022                                        const int run_number,
1023                                        const int event_number)
1024 {
1025   ++m_counter_raw_pairs;
1026
1027   if (!m_write_same_sign_pairs && track1.charge == track2.charge)
1028   {
1029     ++m_counter_reject_charge;
1030     return false;
1031   }
1032
1033   if (!passes_preselection(track1, track2, primary_vertex))
1034   {
1035     ++m_counter_reject_preselection;
1036     return false;
1037   }
1038
1039   Vec3 pca1;
1040   Vec3 pca2;
1041   Vec3 mom1 = track1.momentum;
1042   Vec3 mom2 = track2.momentum;
1043   double pair_dca = 0.0;
1044   double theta1 = quiet_nan();
1045   double theta2 = quiet_nan();
1046   std::pair<double, double> dca1;
1047   std::pair<double, double> dca2;
1048
1049   if (m_fit_kalman_tracks && track1.has_kalman && track2.has_kalman)
1050   {
1051     const auto pca_start = std::chrono::steady_clock::now();
1052     auto candidates = kalman_pca_candidates(
1053         track1.kalman, track2.kalman, m_kalman_config, primary_vertex,
1054         m_kalman_max_upstream_cm, m_kalman_downstream_margin_cm,
1055         m_coarse_steps, m_pca_candidates);
1056     m_timing_kalman_pca_seconds += std::chrono::duration<double>(
1057                                        std::chrono::steady_clock::now() - pca_start)
1058                                        .count();
1059     if (candidates.empty())
1060     {
1061       ++m_counter_reject_pca;
1062       return false;
1063     }
1064
1065     auto best = candidates.front();
1066     if (m_prefer_positive_pointing)
1067     {
1068       double best_score = std::numeric_limits<double>::max();
1069       for (const auto &candidate : candidates)
1070       {
1071         const Vec3 cand_mom1 = kalman_momentum(track1.kalman, candidate.s1, m_kalman_config, primary_vertex);
1072         const Vec3 cand_mom2 = kalman_momentum(track2.kalman, candidate.s2, m_kalman_config, primary_vertex);
1073         const Vec3 cand_vertex = scale(add(candidate.pca1, candidate.pca2), 0.5);
1074         const Vec3 flight = subtract(cand_vertex, primary_vertex);
1075         const Vec3 total_mom = add(cand_mom1, cand_mom2);
1076         const double cos_theta = vector_cosine(flight, total_mom);
1077         const double penalty = (std::isfinite(cos_theta) && cos_theta > 0.0) ? 0.0 : 1000.0;
1078         const double score = penalty + candidate.dca - 1e-3 * cos_theta;
1079         if (score < best_score)
1080         {
1081           best_score = score;
1082           best = candidate;
1083         }
1084       }
1085     }
1086
1087     pca1 = best.pca1;
1088     pca2 = best.pca2;
1089     pair_dca = best.dca;
1090     theta1 = best.s1;
1091     theta2 = best.s2;
1092     mom1 = kalman_momentum(track1.kalman, theta1, m_kalman_config, primary_vertex);
1093     mom2 = kalman_momentum(track2.kalman, theta2, m_kalman_config, primary_vertex);
1094     dca1 = track1.has_vertex_dca
1095                ? track1.vertex_dca
1096                : TpcTrackKalmanFitter::dca_to_vertex(
1097                      track1.kalman, primary_vertex, &m_kalman_config);
1098     dca2 = track2.has_vertex_dca
1099                ? track2.vertex_dca
1100                : TpcTrackKalmanFitter::dca_to_vertex(
1101                      track2.kalman, primary_vertex, &m_kalman_config);
1102   }
1103   else if (m_fit_helix_tracks && track1.has_helix && track2.has_helix)
1104   {
1105     std::vector<HelixPca> candidates;
1106     if (track1.has_helix_search_range && track2.has_helix_search_range)
1107     {
1108       candidates = TpcTrackHelixFitter::pca_candidates_in_ranges(
1109           track1.helix, track2.helix,
1110           track1.helix_search_range, track2.helix_search_range,
1111           m_coarse_steps, m_pca_candidates);
1112     }
1113     else
1114     {
1115       candidates = helix_helix_pca_candidates(
1116           track1.helix, track2.helix, m_theta_extension, m_coarse_steps,
1117           m_downstream_margin, m_pca_candidates);
1118     }
1119     if (candidates.empty())
1120     {
1121       ++m_counter_reject_pca;
1122       return false;
1123     }
1124
1125     auto best = candidates.front();
1126     if (m_prefer_positive_pointing)
1127     {
1128       double best_score = std::numeric_limits<double>::max();
1129       for (const auto &candidate : candidates)
1130       {
1131         const Vec3 cand_mom1 = helix_momentum(track1.helix, candidate.theta1);
1132         const Vec3 cand_mom2 = helix_momentum(track2.helix, candidate.theta2);
1133         const Vec3 cand_vertex = scale(add(candidate.pca1, candidate.pca2), 0.5);
1134         const Vec3 flight = subtract(cand_vertex, primary_vertex);
1135         const Vec3 total_mom = add(cand_mom1, cand_mom2);
1136         const double cos_theta = vector_cosine(flight, total_mom);
1137         const double penalty = (std::isfinite(cos_theta) && cos_theta > 0.0) ? 0.0 : 1000.0;
1138         const double score = penalty + candidate.dca - 1e-3 * cos_theta;
1139         if (score < best_score)
1140         {
1141           best_score = score;
1142           best = candidate;
1143         }
1144       }
1145     }
1146
1147     pca1 = best.pca1;
1148     pca2 = best.pca2;
1149     pair_dca = best.dca;
1150     theta1 = best.theta1;
1151     theta2 = best.theta2;
1152     mom1 = helix_momentum(track1.helix, theta1);
1153     mom2 = helix_momentum(track2.helix, theta2);
1154     dca1 = helix_dca_to_vertex(track1.helix, primary_vertex);
1155     dca2 = helix_dca_to_vertex(track2.helix, primary_vertex);
1156   }
1157   else
1158   {
1159     LinePca pca;
1160     if (!line_line_pca(track1.position, track1.momentum, track2.position, track2.momentum, pca, true))
1161     {
1162       ++m_counter_reject_pca;
1163       return false;
1164     }
1165     pca1 = pca.pca1;
1166     pca2 = pca.pca2;
1167     pair_dca = pca.dca;
1168     dca1 = track_dca_to_vertex(track1.position, track1.momentum, primary_vertex);
1169     dca2 = track_dca_to_vertex(track2.position, track2.momentum, primary_vertex);
1170   }
1171
1172   const Vec3 pair_vertex = scale(add(pca1, pca2), 0.5);
1173   const Vec3 total_mom = add(mom1, mom2);
1174   const Vec3 flight = subtract(pair_vertex, primary_vertex);
1175   const double cos_theta = vector_cosine(flight, total_mom);
1176   if (!std::isfinite(cos_theta) || norm(flight) <= 0.0 || norm(total_mom) <= 0.0)
1177   {
1178     ++m_counter_reject_pointing;
1179     return false;
1180   }
1181
1182   const Vec3 &pplus = (track1.charge > 0) ? mom1 : mom2;
1183   const Vec3 &pminus = (track1.charge > 0) ? mom2 : mom1;
1184   double alpha = 0.0;
1185   double qt = 0.0;
1186   if (!armenteros(pplus, pminus, alpha, qt))
1187   {
1188     ++m_counter_reject_ap;
1189     return false;
1190   }
1191
1192   if (!passes_pair_selection(pca1, pca2, pair_vertex, primary_vertex,
1193                              pair_dca, cos_theta, alpha))
1194   {
1195     ++m_counter_reject_pair_selection;
1196     return false;
1197   }
1198
1199   const Vec3 &truth_pplus = (track1.charge > 0) ? track1.truth_momentum : track2.truth_momentum;
1200   const Vec3 &truth_pminus = (track1.charge > 0) ? track2.truth_momentum : track1.truth_momentum;
1201   double truth_alpha = quiet_nan();
1202   double truth_qt = quiet_nan();
1203   armenteros(truth_pplus, truth_pminus, truth_alpha, truth_qt);
1204
1205   Vec3 true_decay{quiet_nan(), quiet_nan(), quiet_nan()};
1206   double pca_to_true_3d = quiet_nan();
1207   double pca_to_true_xy = quiet_nan();
1208   double pca_to_true_z = quiet_nan();
1209   if (track1.parent_id != 0 && track1.parent_id == track2.parent_id)
1210   {
1211     true_decay = scale(add(track1.truth_vertex, track2.truth_vertex), 0.5);
1212     const Vec3 delta = subtract(pair_vertex, true_decay);
1213     pca_to_true_3d = norm(delta);
1214     pca_to_true_xy = std::sqrt(square(delta.x) + square(delta.y));
1215     pca_to_true_z = std::abs(delta.z);
1216   }
1217
1218   const Vec3 positive_mom = (track1.charge > 0) ? mom1 : mom2;
1219   const Vec3 negative_mom = (track1.charge > 0) ? mom2 : mom1;
1220
1221   reset_pair_row();
1222   m_pair.run = run_number;
1223   m_pair.evt = event_number;
1224   m_pair.cross1 = 0;
1225   m_pair.cross2 = 0;
1226   m_pair.px1 = static_cast<float>(mom1.x);
1227   m_pair.py1 = static_cast<float>(mom1.y);
1228   m_pair.pz1 = static_cast<float>(mom1.z);
1229   m_pair.px2 = static_cast<float>(mom2.x);
1230   m_pair.py2 = static_cast<float>(mom2.y);
1231   m_pair.pz2 = static_cast<float>(mom2.z);
1232   m_pair.dca_xy1 = static_cast<float>(dca1.first);
1233   m_pair.dca_z1 = static_cast<float>(dca1.second);
1234   m_pair.dca_xy2 = static_cast<float>(dca2.first);
1235   m_pair.dca_z2 = static_cast<float>(dca2.second);
1236   m_pair.pairDCA = static_cast<float>(pair_dca);
1237   m_pair.alpha = static_cast<float>(alpha);
1238   m_pair.qT = static_cast<float>(qt);
1239   m_pair.charge1 = static_cast<float>(track1.charge);
1240   m_pair.charge2 = static_cast<float>(track2.charge);
1241   m_pair.dedx_1 = track1.has_dedx ? static_cast<float>(track1.dedx) : quiet_nan();
1242   m_pair.dedx_2 = track2.has_dedx ? static_cast<float>(track2.dedx) : quiet_nan();
1243   m_pair.cosThetaReco = static_cast<float>(cos_theta);
1244   m_pair.Lproj = static_cast<float>(norm(flight));
1245
1246   m_pair.pca_x = static_cast<float>(pair_vertex.x);
1247   m_pair.pca_y = static_cast<float>(pair_vertex.y);
1248   m_pair.pca_z = static_cast<float>(pair_vertex.z);
1249   m_pair.pca1_x = static_cast<float>(pca1.x);
1250   m_pair.pca1_y = static_cast<float>(pca1.y);
1251   m_pair.pca1_z = static_cast<float>(pca1.z);
1252   m_pair.pca2_x = static_cast<float>(pca2.x);
1253   m_pair.pca2_y = static_cast<float>(pca2.y);
1254   m_pair.pca2_z = static_cast<float>(pca2.z);
1255
1256   m_pair.v0_px = static_cast<float>(total_mom.x);
1257   m_pair.v0_py = static_cast<float>(total_mom.y);
1258   m_pair.v0_pz = static_cast<float>(total_mom.z);
1259   m_pair.v0_pt = static_cast<float>(pt(total_mom));
1260   m_pair.mass_Kshort = static_cast<float>(invariant_mass(mom1, kPionMass, mom2, kPionMass));
1261   m_pair.mass_Lambda = static_cast<float>(invariant_mass(positive_mom, kProtonMass, negative_mom, kPionMass));
1262   m_pair.mass_AntiLambda = static_cast<float>(invariant_mass(positive_mom, kPionMass, negative_mom, kProtonMass));
1263
1264   m_pair.true_decay_x = static_cast<float>(true_decay.x);
1265   m_pair.true_decay_y = static_cast<float>(true_decay.y);
1266   m_pair.true_decay_z = static_cast<float>(true_decay.z);
1267   m_pair.pca_to_true_3d = static_cast<float>(pca_to_true_3d);
1268   m_pair.pca_to_true_xy = static_cast<float>(pca_to_true_xy);
1269   m_pair.pca_to_true_z = static_cast<float>(pca_to_true_z);
1270   m_pair.truth_alpha = static_cast<float>(truth_alpha);
1271   m_pair.truth_qT = static_cast<float>(truth_qt);
1272   m_pair.delta_alpha = static_cast<float>(alpha - truth_alpha);
1273   m_pair.delta_qT = static_cast<float>(qt - truth_qt);
1274   m_pair.truth_px1 = static_cast<float>(track1.truth_momentum.x);
1275   m_pair.truth_py1 = static_cast<float>(track1.truth_momentum.y);
1276   m_pair.truth_pz1 = static_cast<float>(track1.truth_momentum.z);
1277   m_pair.truth_px2 = static_cast<float>(track2.truth_momentum.x);
1278   m_pair.truth_py2 = static_cast<float>(track2.truth_momentum.y);
1279   m_pair.truth_pz2 = static_cast<float>(track2.truth_momentum.z);
1280   m_pair.cos_mom1_truth = static_cast<float>(vector_cosine(mom1, track1.truth_momentum));
1281   m_pair.cos_mom2_truth = static_cast<float>(vector_cosine(mom2, track2.truth_momentum));
1282   m_pair.pca_theta1 = static_cast<float>(theta1);
1283   m_pair.pca_theta2 = static_cast<float>(theta2);
1284   if (track1.has_kalman)
1285   {
1286     m_pair.kalman_chi2_1 = static_cast<float>(track1.kalman.chi2);
1287     m_pair.kalman_ndof1 = track1.kalman.ndof;
1288     m_pair.kalman_chi2_ndf1 =
1289         (track1.kalman.ndof > 0 && std::isfinite(track1.kalman.chi2))
1290             ? static_cast<float>(track1.kalman.chi2 / track1.kalman.ndof)
1291             : quiet_nan();
1292   }
1293   if (track2.has_kalman)
1294   {
1295     m_pair.kalman_chi2_2 = static_cast<float>(track2.kalman.chi2);
1296     m_pair.kalman_ndof2 = track2.kalman.ndof;
1297     m_pair.kalman_chi2_ndf2 =
1298         (track2.kalman.ndof > 0 && std::isfinite(track2.kalman.chi2))
1299             ? static_cast<float>(track2.kalman.chi2 / track2.kalman.ndof)
1300             : quiet_nan();
1301   }
1302   m_pair.quality1 = std::isfinite(track1.fit_chi2_ndf)
1303                         ? static_cast<float>(track1.fit_chi2_ndf)
1304                         : quiet_nan();
1305   m_pair.quality2 = std::isfinite(track2.fit_chi2_ndf)
1306                         ? static_cast<float>(track2.fit_chi2_ndf)
1307                         : quiet_nan();
1308   m_pair.track_id1 = track1.track_id;
1309   m_pair.track_id2 = track2.track_id;
1310   m_pair.pid1 = track1.pid;
1311   m_pair.pid2 = track2.pid;
1312   m_pair.parent_id1 = track1.parent_id;
1313   m_pair.parent_id2 = track2.parent_id;
1314   m_pair.parent_pid = (track1.parent_id != 0 && track1.parent_id == track2.parent_id) ? track1.parent_pid : 0;
1315   m_pair.npoints1 = static_cast<short>(track1.npoints);
1316   m_pair.npoints2 = static_cast<short>(track2.npoints);
1317
1318   m_pair_tree->Fill();
1319   ++m_counter_written;
1320   return true;
1321 }
1322
1323 void TpcV0CandidateTree::fill_track_row(const Tracklet &tracklet,
1324                                         const Vec3 &primary_vertex,
1325                                         const int run_number,
1326                                         const int event_number)
1327 {
1328   reset_track_row();
1329   const float nan = quiet_nan();
1330
1331   m_track.run = run_number;
1332   m_track.evt = event_number;
1333   m_track.track_id = tracklet.track_id;
1334   m_track.shower_id = tracklet.shower_id;
1335   m_track.pid = tracklet.pid;
1336   m_track.parent_id = tracklet.parent_id;
1337   m_track.parent_pid = tracklet.parent_pid;
1338   m_track.charge = static_cast<double>(tracklet.charge);
1339   m_track.side = tracklet.side;
1340   m_track.npoints = tracklet.npoints;
1341   m_track.ntpc_clusters = tracklet.ntpc_clusters;
1342   m_track.has_helix = tracklet.has_helix ? 1 : 0;
1343   m_track.has_kalman = tracklet.has_kalman ? 1 : 0;
1344   m_track.is_primary = tracklet.is_primary;
1345
1346   if (m_kalman_config.collect_innovation_components)
1347   {
1348     m_track.kalman_measurement_sigma_r =
1349         static_cast<float>(m_kalman_config.meas_sigma_r_cm);
1350     m_track.kalman_measurement_sigma_rphi =
1351         static_cast<float>(m_kalman_config.meas_sigma_rphi_cm);
1352     m_track.kalman_measurement_sigma_z =
1353         static_cast<float>(m_kalman_config.meas_sigma_z_cm);
1354   }
1355
1356   Vec3 row_position = tracklet.position;
1357   Vec3 row_momentum = tracklet.momentum;
1358   std::array<double, TpcTrackKalmanFitter::StateDim> kalman_row_state{};
1359   bool has_kalman_row_state = false;
1360   if (tracklet.has_kalman)
1361   {
1362     kalman_row_state = TpcTrackKalmanFitter::propagation_state(tracklet.kalman, primary_vertex);
1363     row_position = TpcTrackKalmanFitter::state_position(kalman_row_state);
1364     row_momentum = TpcTrackKalmanFitter::state_momentum(kalman_row_state);
1365     has_kalman_row_state = true;
1366   }
1367
1368   m_track.px = row_momentum.x;
1369   m_track.py = row_momentum.y;
1370   m_track.pz = row_momentum.z;
1371   m_track.pt = pt(row_momentum);
1372   m_track.p = norm(row_momentum);
1373   m_track.eta = m_track.pt > 0.0
1374                     ? std::asinh(row_momentum.z / m_track.pt)
1375                     : static_cast<double>(nan);
1376   m_track.dedx = tracklet.has_dedx ? tracklet.dedx : static_cast<double>(nan);
1377   m_track.x = static_cast<float>(row_position.x);
1378   m_track.y = static_cast<float>(row_position.y);
1379   m_track.z = static_cast<float>(row_position.z);
1380
1381   if (!tracklet.points.empty())
1382   {
1383     const auto &first = tracklet.points.front().position;
1384     const auto &last = tracklet.points.back().position;
1385     m_track.first_x = static_cast<float>(first.x);
1386     m_track.first_y = static_cast<float>(first.y);
1387     m_track.first_z = static_cast<float>(first.z);
1388     m_track.first_r = static_cast<float>(pt(first));
1389     m_track.last_x = static_cast<float>(last.x);
1390     m_track.last_y = static_cast<float>(last.y);
1391     m_track.last_z = static_cast<float>(last.z);
1392     m_track.last_r = static_cast<float>(pt(last));
1393   }
1394
1395   m_track.cluster_index.reserve(tracklet.points.size());
1396   m_track.cluster_side.reserve(tracklet.points.size());
1397   m_track.layer.reserve(tracklet.points.size());
1398   m_track.cluster_z.reserve(tracklet.points.size());
1399   m_track.cluster_r.reserve(tracklet.points.size());
1400   m_track.cluster_phi.reserve(tracklet.points.size());
1401   m_track.residual_z.reserve(tracklet.points.size());
1402   m_track.residual_r.reserve(tracklet.points.size());
1403   m_track.residual_rphi.reserve(tracklet.points.size());
1404
1405   double previous_theta = tracklet.has_helix ? tracklet.helix.theta_first : 0.0;
1406   bool have_previous_theta = false;
1407   for (std::size_t index = 0; index < tracklet.points.size(); ++index)
1408   {
1409     const auto &point = tracklet.points[index];
1410     const Vec3 &cluster = point.position;
1411     const double cluster_r = pt(cluster);
1412     const double cluster_phi = std::atan2(cluster.y, cluster.x);
1413
1414     double residual_r = std::numeric_limits<double>::quiet_NaN();
1415     double residual_rphi = std::numeric_limits<double>::quiet_NaN();
1416     double residual_z = std::numeric_limits<double>::quiet_NaN();
1417
1418     if (tracklet.has_kalman && index < tracklet.kalman.states_smoothed.size())
1419     {
1420       const Vec3 fit_position = TpcTrackKalmanFitter::state_position(tracklet.kalman.states_smoothed[index]);
1421       const Vec3 delta = subtract(cluster, fit_position);
1422       const double fit_phi = std::atan2(fit_position.y, fit_position.x);
1423       residual_r = std::cos(fit_phi) * delta.x + std::sin(fit_phi) * delta.y;
1424       residual_rphi = -std::sin(fit_phi) * delta.x + std::cos(fit_phi) * delta.y;
1425       residual_z = delta.z;
1426     }
1427     else if (tracklet.has_helix)
1428     {
1429       double theta_hint = std::atan2(cluster.y - tracklet.helix.cy,
1430                                      cluster.x - tracklet.helix.cx);
1431       theta_hint = unwrap_to_near(theta_hint,
1432                                   have_previous_theta ? previous_theta
1433                                                       : tracklet.helix.theta_first);
1434
1435       double theta = theta_hint;
1436       Vec3 fit_position;
1437       if (!helix_point_at_beam_radius(tracklet.helix, cluster_r, theta_hint, theta, fit_position))
1438       {
1439         fit_position = helix_point(tracklet.helix, theta_hint);
1440         theta = theta_hint;
1441       }
1442
1443       if (finite(fit_position))
1444       {
1445         previous_theta = theta;
1446         have_previous_theta = true;
1447         const double fit_r = pt(fit_position);
1448         const double fit_phi = std::atan2(fit_position.y, fit_position.x);
1449         residual_r = cluster_r - fit_r;
1450         residual_rphi = cluster_r * normalize_phi(cluster_phi - fit_phi);
1451         residual_z = cluster.z - fit_position.z;
1452       }
1453     }
1454
1455     m_track.cluster_index.push_back(static_cast<unsigned int>(index));
1456     m_track.cluster_side.push_back(tracklet.side);
1457     m_track.layer.push_back(static_cast<unsigned int>(std::max(point.layer, 0)));
1458     m_track.cluster_z.push_back(cluster.z);
1459     m_track.cluster_r.push_back(cluster_r);
1460     m_track.cluster_phi.push_back(cluster_phi);
1461     m_track.residual_z.push_back(std::isfinite(residual_z) ? residual_z : static_cast<double>(nan));
1462     m_track.residual_r.push_back(std::isfinite(residual_r) ? residual_r : static_cast<double>(nan));
1463     m_track.residual_rphi.push_back(std::isfinite(residual_rphi) ? residual_rphi : static_cast<double>(nan));
1464   }
1465
1466   const Vec3 &track_vertex = tracklet.has_pattern_vertex
1467                                  ? tracklet.pattern_vertex
1468                                  : primary_vertex;
1469   auto dca = tracklet.vertex_dca;
1470   if (!tracklet.has_vertex_dca)
1471   {
1472     dca = fitted_track_dca_to_vertex(tracklet, track_vertex);
1473   }
1474   m_track.dca_xy = static_cast<float>(dca.first);
1475   m_track.dca_z = static_cast<float>(dca.second);
1476   m_track.vertex_x = track_vertex.x;
1477   m_track.vertex_y = track_vertex.y;
1478   m_track.vertex_z = track_vertex.z;
1479   m_track.vertex_from_upstream = tracklet.has_pattern_vertex ? 1 : 0;
1480   if (tracklet.has_pattern_vertex)
1481   {
1482     m_track.vertex_z_rms = tracklet.pattern_vertex_z_rms;
1483     m_track.vertex_ntracks = tracklet.pattern_vertex_ntracks;
1484   }
1485   if (tracklet.has_beamline_pca)
1486   {
1487     m_track.pca_x = tracklet.beamline_pca.x;
1488     m_track.pca_y = tracklet.beamline_pca.y;
1489     m_track.pca_z = tracklet.beamline_pca.z;
1490     m_track.rDCA_zero = tracklet.rdca_zero;
1491     m_track.zDCA = tracklet.beamline_pca.z - track_vertex.z;
1492   }
1493
1494   if (tracklet.has_helix)
1495   {
1496     m_track.helix_cx = static_cast<float>(tracklet.helix.cx);
1497     m_track.helix_cy = static_cast<float>(tracklet.helix.cy);
1498     m_track.helix_radius = static_cast<float>(tracklet.helix.radius);
1499     m_track.helix_z0 = static_cast<float>(tracklet.helix.z0);
1500     m_track.helix_pitch = static_cast<float>(tracklet.helix.pitch);
1501     m_track.helix_theta_first = static_cast<float>(tracklet.helix.theta_first);
1502     m_track.helix_theta_last = static_cast<float>(tracklet.helix.theta_last);
1503     m_track.helix_direction = static_cast<float>(tracklet.helix.direction);
1504     if (tracklet.has_helix_search_range)
1505     {
1506       const auto &range = tracklet.helix_search_range;
1507       m_track.helix_search_anchored = 1;
1508       m_track.helix_anchor_point_index = range.anchor_point_index;
1509       m_track.helix_anchor_theta = static_cast<float>(range.anchor_theta);
1510       m_track.helix_anchor_path_cm = static_cast<float>(range.anchor_path_cm);
1511       m_track.helix_anchor_residual_cm = static_cast<float>(range.anchor_residual_cm);
1512       m_track.helix_search_theta_min = static_cast<float>(range.theta_min);
1513       m_track.helix_search_theta_max = static_cast<float>(range.theta_max);
1514       m_track.helix_search_upstream_cm = static_cast<float>(range.upstream_cm);
1515       m_track.helix_search_downstream_cm = static_cast<float>(range.downstream_cm);
1516     }
1517   }
1518
1519   if (tracklet.has_kalman)
1520   {
1521     m_track.kalman_chi2 = static_cast<float>(tracklet.kalman.chi2);
1522     m_track.kalman_ndof = tracklet.kalman.ndof;
1523     m_track.kalman_naccepted = static_cast<unsigned int>(tracklet.kalman.naccepted);
1524     m_track.kalman_nrejected = static_cast<unsigned int>(tracklet.kalman.nrejected);
1525     m_track.kalman_measurement_chi2 = tracklet.kalman.measurement_chi2;
1526     m_track.kalman_measurement_used = tracklet.kalman.measurement_used;
1527     if (m_kalman_config.collect_innovation_components)
1528     {
1529       m_track.kalman_measurement_in_seed = tracklet.kalman.measurement_in_seed;
1530       m_track.kalman_innovation_residual_r = tracklet.kalman.innovation_residual_r;
1531       m_track.kalman_innovation_residual_rphi = tracklet.kalman.innovation_residual_rphi;
1532       m_track.kalman_innovation_residual_z = tracklet.kalman.innovation_residual_z;
1533       m_track.kalman_prediction_sigma_r = tracklet.kalman.prediction_sigma_r;
1534       m_track.kalman_prediction_sigma_rphi = tracklet.kalman.prediction_sigma_rphi;
1535       m_track.kalman_prediction_sigma_z = tracklet.kalman.prediction_sigma_z;
1536       m_track.kalman_innovation_sigma_r = tracklet.kalman.innovation_sigma_r;
1537       m_track.kalman_innovation_sigma_rphi = tracklet.kalman.innovation_sigma_rphi;
1538       m_track.kalman_innovation_sigma_z = tracklet.kalman.innovation_sigma_z;
1539       m_track.kalman_innovation_rho_r_rphi = tracklet.kalman.innovation_rho_r_rphi;
1540       m_track.kalman_innovation_rho_r_z = tracklet.kalman.innovation_rho_r_z;
1541       m_track.kalman_innovation_rho_rphi_z = tracklet.kalman.innovation_rho_rphi_z;
1542       m_track.kalman_innovation_whitened_0 = tracklet.kalman.innovation_whitened_0;
1543       m_track.kalman_innovation_whitened_1 = tracklet.kalman.innovation_whitened_1;
1544       m_track.kalman_innovation_whitened_2 = tracklet.kalman.innovation_whitened_2;
1545     }
1546     if (has_kalman_row_state)
1547     {
1548       const auto &state = kalman_row_state;
1549       const double qop_t = state[TpcTrackKalmanFitter::QOverPt];
1550       const double omega = 0.003 * tracklet.kalman.bfield_t * qop_t;
1551       m_track.kalman_qop_t = static_cast<float>(qop_t);
1552       m_track.kalman_omega = static_cast<float>(omega);
1553       if (std::abs(omega) > 1.0e-12)
1554       {
1555         const double center_x = state[TpcTrackKalmanFitter::X] -
1556                                 std::sin(state[TpcTrackKalmanFitter::Phi]) / omega;
1557         const double center_y = state[TpcTrackKalmanFitter::Y] +
1558                                 std::cos(state[TpcTrackKalmanFitter::Phi]) / omega;
1559         m_track.kalman_cx = static_cast<float>(center_x);
1560         m_track.kalman_cy = static_cast<float>(center_y);
1561         m_track.kalman_radius = static_cast<float>(std::abs(1.0 / omega));
1562       }
1563     }
1564   }
1565
1566   m_track.fit_chi2 = std::isfinite(tracklet.fit_chi2)
1567                          ? static_cast<float>(tracklet.fit_chi2)
1568                          : nan;
1569   m_track.fit_ndf = tracklet.fit_ndf;
1570   m_track.quality = std::isfinite(tracklet.fit_chi2_ndf)
1571                         ? static_cast<float>(tracklet.fit_chi2_ndf)
1572                         : nan;
1573
1574   m_track.truth_px = static_cast<float>(tracklet.truth_momentum.x);
1575   m_track.truth_py = static_cast<float>(tracklet.truth_momentum.y);
1576   m_track.truth_pz = static_cast<float>(tracklet.truth_momentum.z);
1577   m_track.cos_mom_truth = static_cast<float>(vector_cosine(row_momentum, tracklet.truth_momentum));
1578
1579   m_track_tree->Fill();
1580   ++m_counter_tracks_written;
1581 }
1582
1583 void TpcV0CandidateTree::fill_cluster_residual_rows(const Tracklet &tracklet,
1584                                                     const Vec3 & /*primary_vertex*/,
1585                                                     const int run_number,
1586                                                     const int event_number)
1587 {
1588   if (!m_cluster_residual_tree || tracklet.points.empty())
1589   {
1590     return;
1591   }
1592
1593   const float nan = quiet_nan();
1594   double fit_chi2 = std::numeric_limits<double>::quiet_NaN();
1595   int fit_ndf = 0;
1596
1597   const double sigma_rphi = std::max(m_kalman_config.meas_sigma_rphi_cm,
1598                                      m_kalman_config.min_measurement_sigma_cm);
1599   const double sigma_z = std::max(m_kalman_config.meas_sigma_z_cm,
1600                                   m_kalman_config.min_measurement_sigma_cm);
1601
1602   auto helix_residual = [&](const TruthPoint &point,
1603                             double &previous_theta,
1604                             bool &have_previous_theta,
1605                             Vec3 &fit_position,
1606                             double &residual_r,
1607                             double &residual_rphi,
1608                             double &residual_z) -> bool
1609   {
1610     if (!tracklet.has_helix)
1611     {
1612       return false;
1613     }
1614
1615     const Vec3 &cluster = point.position;
1616     const double cluster_r = pt(cluster);
1617     double theta_hint = std::atan2(cluster.y - tracklet.helix.cy,
1618                                    cluster.x - tracklet.helix.cx);
1619     theta_hint = unwrap_to_near(theta_hint,
1620                                 have_previous_theta ? previous_theta
1621                                                     : tracklet.helix.theta_first);
1622
1623     double theta = theta_hint;
1624     if (!helix_point_at_beam_radius(tracklet.helix, cluster_r, theta_hint, theta, fit_position))
1625     {
1626       fit_position = helix_point(tracklet.helix, theta_hint);
1627       theta = theta_hint;
1628       if (!finite(fit_position))
1629       {
1630         return false;
1631       }
1632     }
1633
1634     previous_theta = theta;
1635     have_previous_theta = true;
1636
1637     const double fit_r = pt(fit_position);
1638     const double cluster_phi = std::atan2(cluster.y, cluster.x);
1639     const double fit_phi = std::atan2(fit_position.y, fit_position.x);
1640     residual_r = cluster_r - fit_r;
1641     residual_rphi = cluster_r * normalize_phi(cluster_phi - fit_phi);
1642     residual_z = cluster.z - fit_position.z;
1643     return std::isfinite(residual_r) &&
1644            std::isfinite(residual_rphi) &&
1645            std::isfinite(residual_z);
1646   };
1647
1648   if (tracklet.has_kalman)
1649   {
1650     fit_chi2 = tracklet.kalman.chi2;
1651     fit_ndf = tracklet.kalman.ndof;
1652   }
1653   else if (tracklet.has_helix)
1654   {
1655     double chi2 = 0.0;
1656     int nresiduals = 0;
1657     double previous_theta = tracklet.helix.theta_first;
1658     bool have_previous_theta = false;
1659     for (const auto &point : tracklet.points)
1660     {
1661       Vec3 fit_position;
1662       double residual_r = 0.0;
1663       double residual_rphi = 0.0;
1664       double residual_z = 0.0;
1665       if (!helix_residual(point, previous_theta, have_previous_theta,
1666                           fit_position, residual_r, residual_rphi, residual_z))
1667       {
1668         continue;
1669       }
1670       chi2 += square(residual_rphi / sigma_rphi) + square(residual_z / sigma_z);
1671       nresiduals += 2;
1672     }
1673     fit_chi2 = chi2;
1674     fit_ndf = std::max(0, nresiduals - 5);
1675   }
1676
1677   double previous_theta = tracklet.has_helix ? tracklet.helix.theta_first : 0.0;
1678   bool have_previous_theta = false;
1679   for (std::size_t index = 0; index < tracklet.points.size(); ++index)
1680   {
1681     const auto &point = tracklet.points[index];
1682     const Vec3 &cluster = point.position;
1683     Vec3 fit_position{nan, nan, nan};
1684     double residual_r = std::numeric_limits<double>::quiet_NaN();
1685     double residual_rphi = std::numeric_limits<double>::quiet_NaN();
1686     double residual_z = std::numeric_limits<double>::quiet_NaN();
1687
1688     if (tracklet.has_kalman && index < tracklet.kalman.states_smoothed.size())
1689     {
1690       fit_position = TpcTrackKalmanFitter::state_position(tracklet.kalman.states_smoothed[index]);
1691       const Vec3 delta = subtract(cluster, fit_position);
1692       const double fit_phi = std::atan2(fit_position.y, fit_position.x);
1693       residual_r = std::cos(fit_phi) * delta.x + std::sin(fit_phi) * delta.y;
1694       residual_rphi = -std::sin(fit_phi) * delta.x + std::cos(fit_phi) * delta.y;
1695       residual_z = delta.z;
1696     }
1697     else if (tracklet.has_helix)
1698     {
1699       helix_residual(point, previous_theta, have_previous_theta,
1700                      fit_position, residual_r, residual_rphi, residual_z);
1701     }
1702
1703     reset_cluster_residual_row();
1704     m_cluster_residual.run = run_number;
1705     m_cluster_residual.evt = event_number;
1706     m_cluster_residual.track_id = tracklet.track_id;
1707     m_cluster_residual.charge = tracklet.charge;
1708     m_cluster_residual.side = (cluster.z >= 0.0) ? 1 : -1;
1709     m_cluster_residual.layer = point.layer;
1710     m_cluster_residual.cluster_index = static_cast<int>(index);
1711     m_cluster_residual.ntp_cluster = tracklet.npoints;
1712     m_cluster_residual.npoints = tracklet.npoints;
1713     m_cluster_residual.has_helix = tracklet.has_helix ? 1 : 0;
1714     m_cluster_residual.has_kalman = tracklet.has_kalman ? 1 : 0;
1715
1716     m_cluster_residual.cluster_x = static_cast<float>(cluster.x);
1717     m_cluster_residual.cluster_y = static_cast<float>(cluster.y);
1718     m_cluster_residual.cluster_z = static_cast<float>(cluster.z);
1719     m_cluster_residual.cluster_r = static_cast<float>(pt(cluster));
1720     m_cluster_residual.cluster_phi = static_cast<float>(std::atan2(cluster.y, cluster.x));
1721
1722     m_cluster_residual.fit_x = static_cast<float>(fit_position.x);
1723     m_cluster_residual.fit_y = static_cast<float>(fit_position.y);
1724     m_cluster_residual.fit_z = static_cast<float>(fit_position.z);
1725     m_cluster_residual.fit_r = static_cast<float>(pt(fit_position));
1726     m_cluster_residual.fit_phi = static_cast<float>(std::atan2(fit_position.y, fit_position.x));
1727
1728     m_cluster_residual.residual_x = static_cast<float>(cluster.x - fit_position.x);
1729     m_cluster_residual.residual_y = static_cast<float>(cluster.y - fit_position.y);
1730     m_cluster_residual.residual_z = static_cast<float>(residual_z);
1731     m_cluster_residual.residual_r = static_cast<float>(residual_r);
1732     m_cluster_residual.residual_rphi = static_cast<float>(residual_rphi);
1733
1734     m_cluster_residual.fit_chi2 = static_cast<float>(fit_chi2);
1735     m_cluster_residual.fit_ndf = fit_ndf;
1736     m_cluster_residual.fit_chi2_ndf =
1737         (fit_ndf > 0 && std::isfinite(fit_chi2)) ? static_cast<float>(fit_chi2 / fit_ndf) : nan;
1738
1739     m_cluster_residual_tree->Fill();
1740     ++m_counter_cluster_residuals_written;
1741   }
1742 }
1743
1744 void TpcV0CandidateTree::assign_fit_quality(Tracklet &tracklet) const
1745 {
1746   tracklet.fit_chi2 = std::numeric_limits<double>::quiet_NaN();
1747   tracklet.fit_ndf = 0;
1748   tracklet.fit_chi2_ndf = std::numeric_limits<double>::quiet_NaN();
1749
1750   if (tracklet.has_kalman)
1751   {
1752     tracklet.fit_chi2 = tracklet.kalman.chi2;
1753     tracklet.fit_ndf = tracklet.kalman.ndof;
1754   }
1755   else if (tracklet.has_helix)
1756   {
1757     const double sigma_rphi = std::max(m_kalman_config.meas_sigma_rphi_cm,
1758                                        m_kalman_config.min_measurement_sigma_cm);
1759     const double sigma_z = std::max(m_kalman_config.meas_sigma_z_cm,
1760                                     m_kalman_config.min_measurement_sigma_cm);
1761
1762     double chi2 = 0.0;
1763     int nresiduals = 0;
1764     double previous_theta = tracklet.helix.theta_first;
1765     bool have_previous_theta = false;
1766     for (const auto &point : tracklet.points)
1767     {
1768       Vec3 fit_position;
1769       double residual_r = 0.0;
1770       double residual_rphi = 0.0;
1771       double residual_z = 0.0;
1772       if (!helix_cluster_residual(tracklet.helix, point, previous_theta, have_previous_theta,
1773                                   fit_position, residual_r, residual_rphi, residual_z))
1774       {
1775         continue;
1776       }
1777       chi2 += square(residual_rphi / sigma_rphi) + square(residual_z / sigma_z);
1778       nresiduals += 2;
1779     }
1780
1781     tracklet.fit_chi2 = chi2;
1782     tracklet.fit_ndf = std::max(0, nresiduals - 5);
1783   }
1784
1785   if (tracklet.fit_ndf > 0 && std::isfinite(tracklet.fit_chi2))
1786   {
1787     tracklet.fit_chi2_ndf = tracklet.fit_chi2 / tracklet.fit_ndf;
1788   }
1789 }
1790
1791 void TpcV0CandidateTree::reset_pair_row()
1792 {
1793   m_pair = {};
1794   const float nan = quiet_nan();
1795   m_pair.dca_xy1 = nan;
1796   m_pair.dca_z1 = nan;
1797   m_pair.dca_xy2 = nan;
1798   m_pair.dca_z2 = nan;
1799   m_pair.pairDCA = nan;
1800   m_pair.alpha = nan;
1801   m_pair.qT = nan;
1802   m_pair.dedx_1 = nan;
1803   m_pair.dedx_2 = nan;
1804   m_pair.cosThetaReco = nan;
1805   m_pair.Lproj = nan;
1806   m_pair.pca_x = nan;
1807   m_pair.pca_y = nan;
1808   m_pair.pca_z = nan;
1809   m_pair.pca1_x = nan;
1810   m_pair.pca1_y = nan;
1811   m_pair.pca1_z = nan;
1812   m_pair.pca2_x = nan;
1813   m_pair.pca2_y = nan;
1814   m_pair.pca2_z = nan;
1815   m_pair.v0_px = nan;
1816   m_pair.v0_py = nan;
1817   m_pair.v0_pz = nan;
1818   m_pair.v0_pt = nan;
1819   m_pair.mass_Kshort = nan;
1820   m_pair.mass_Lambda = nan;
1821   m_pair.mass_AntiLambda = nan;
1822   m_pair.true_decay_x = nan;
1823   m_pair.true_decay_y = nan;
1824   m_pair.true_decay_z = nan;
1825   m_pair.pca_to_true_3d = nan;
1826   m_pair.pca_to_true_xy = nan;
1827   m_pair.pca_to_true_z = nan;
1828   m_pair.truth_alpha = nan;
1829   m_pair.truth_qT = nan;
1830   m_pair.delta_alpha = nan;
1831   m_pair.delta_qT = nan;
1832   m_pair.truth_px1 = nan;
1833   m_pair.truth_py1 = nan;
1834   m_pair.truth_pz1 = nan;
1835   m_pair.truth_px2 = nan;
1836   m_pair.truth_py2 = nan;
1837   m_pair.truth_pz2 = nan;
1838   m_pair.cos_mom1_truth = nan;
1839   m_pair.cos_mom2_truth = nan;
1840   m_pair.pca_theta1 = nan;
1841   m_pair.pca_theta2 = nan;
1842   m_pair.kalman_chi2_1 = nan;
1843   m_pair.kalman_chi2_2 = nan;
1844   m_pair.kalman_chi2_ndf1 = nan;
1845   m_pair.kalman_chi2_ndf2 = nan;
1846   m_pair.quality1 = nan;
1847   m_pair.quality2 = nan;
1848   m_pair.kalman_ndof1 = 0;
1849   m_pair.kalman_ndof2 = 0;
1850 }
1851
1852 void TpcV0CandidateTree::reset_track_row()
1853 {
1854   m_track = {};
1855   const float nan = quiet_nan();
1856   m_track.px = nan;
1857   m_track.py = nan;
1858   m_track.pz = nan;
1859   m_track.pt = nan;
1860   m_track.p = nan;
1861   m_track.eta = nan;
1862   m_track.dedx = nan;
1863   m_track.x = nan;
1864   m_track.y = nan;
1865   m_track.z = nan;
1866   m_track.first_x = nan;
1867   m_track.first_y = nan;
1868   m_track.first_z = nan;
1869   m_track.first_r = nan;
1870   m_track.last_x = nan;
1871   m_track.last_y = nan;
1872   m_track.last_z = nan;
1873   m_track.last_r = nan;
1874   m_track.dca_xy = nan;
1875   m_track.dca_z = nan;
1876   m_track.vertex_x = nan;
1877   m_track.vertex_y = nan;
1878   m_track.vertex_z = nan;
1879   m_track.vertex_z_rms = nan;
1880   m_track.pca_x = nan;
1881   m_track.pca_y = nan;
1882   m_track.pca_z = nan;
1883   m_track.rDCA_zero = nan;
1884   m_track.zDCA = nan;
1885   m_track.helix_cx = nan;
1886   m_track.helix_cy = nan;
1887   m_track.helix_radius = nan;
1888   m_track.helix_z0 = nan;
1889   m_track.helix_pitch = nan;
1890   m_track.helix_theta_first = nan;
1891   m_track.helix_theta_last = nan;
1892   m_track.helix_direction = nan;
1893   m_track.helix_search_anchored = 0;
1894   m_track.helix_anchor_point_index = -1;
1895   m_track.helix_anchor_theta = nan;
1896   m_track.helix_anchor_path_cm = nan;
1897   m_track.helix_anchor_residual_cm = nan;
1898   m_track.helix_search_theta_min = nan;
1899   m_track.helix_search_theta_max = nan;
1900   m_track.helix_search_upstream_cm = nan;
1901   m_track.helix_search_downstream_cm = nan;
1902   m_track.kalman_chi2 = nan;
1903   m_track.kalman_ndof = 0;
1904   m_track.kalman_qop_t = nan;
1905   m_track.kalman_omega = nan;
1906   m_track.kalman_cx = nan;
1907   m_track.kalman_cy = nan;
1908   m_track.kalman_radius = nan;
1909   m_track.fit_chi2 = nan;
1910   m_track.fit_ndf = 0;
1911   m_track.quality = nan;
1912   m_track.truth_px = nan;
1913   m_track.truth_py = nan;
1914   m_track.truth_pz = nan;
1915   m_track.cos_mom_truth = nan;
1916 }
1917
1918 void TpcV0CandidateTree::reset_cluster_residual_row()
1919 {
1920   m_cluster_residual = {};
1921   const float nan = quiet_nan();
1922   m_cluster_residual.cluster_x = nan;
1923   m_cluster_residual.cluster_y = nan;
1924   m_cluster_residual.cluster_z = nan;
1925   m_cluster_residual.cluster_r = nan;
1926   m_cluster_residual.cluster_phi = nan;
1927   m_cluster_residual.fit_x = nan;
1928   m_cluster_residual.fit_y = nan;
1929   m_cluster_residual.fit_z = nan;
1930   m_cluster_residual.fit_r = nan;
1931   m_cluster_residual.fit_phi = nan;
1932   m_cluster_residual.residual_x = nan;
1933   m_cluster_residual.residual_y = nan;
1934   m_cluster_residual.residual_z = nan;
1935   m_cluster_residual.residual_r = nan;
1936   m_cluster_residual.residual_rphi = nan;
1937   m_cluster_residual.fit_chi2 = nan;
1938   m_cluster_residual.fit_ndf = 0;
1939   m_cluster_residual.fit_chi2_ndf = nan;
1940 }
1941
1942 void TpcV0CandidateTree::create_branches()
1943 {
1944   m_pair_tree->Branch("run", &m_pair.run, "run/I");
1945   m_pair_tree->Branch("evt", &m_pair.evt, "evt/I");
1946   m_pair_tree->Branch("cross1", &m_pair.cross1, "cross1/S");
1947   m_pair_tree->Branch("cross2", &m_pair.cross2, "cross2/S");
1948   m_pair_tree->Branch("px1", &m_pair.px1, "px1/F");
1949   m_pair_tree->Branch("py1", &m_pair.py1, "py1/F");
1950   m_pair_tree->Branch("pz1", &m_pair.pz1, "pz1/F");
1951   m_pair_tree->Branch("px2", &m_pair.px2, "px2/F");
1952   m_pair_tree->Branch("py2", &m_pair.py2, "py2/F");
1953   m_pair_tree->Branch("pz2", &m_pair.pz2, "pz2/F");
1954   m_pair_tree->Branch("dca_xy1", &m_pair.dca_xy1, "dca_xy1/F");
1955   m_pair_tree->Branch("dca_z1", &m_pair.dca_z1, "dca_z1/F");
1956   m_pair_tree->Branch("dca_xy2", &m_pair.dca_xy2, "dca_xy2/F");
1957   m_pair_tree->Branch("dca_z2", &m_pair.dca_z2, "dca_z2/F");
1958   m_pair_tree->Branch("pairDCA", &m_pair.pairDCA, "pairDCA/F");
1959   m_pair_tree->Branch("alpha", &m_pair.alpha, "alpha/F");
1960   m_pair_tree->Branch("qT", &m_pair.qT, "qT/F");
1961   m_pair_tree->Branch("charge1", &m_pair.charge1, "charge1/F");
1962   m_pair_tree->Branch("charge2", &m_pair.charge2, "charge2/F");
1963   m_pair_tree->Branch("dedx_1", &m_pair.dedx_1, "dedx_1/F");
1964   m_pair_tree->Branch("dedx_2", &m_pair.dedx_2, "dedx_2/F");
1965   m_pair_tree->Branch("cosThetaReco", &m_pair.cosThetaReco, "cosThetaReco/F");
1966   m_pair_tree->Branch("Lproj", &m_pair.Lproj, "Lproj/F");
1967   m_pair_tree->Branch("pca_x", &m_pair.pca_x, "pca_x/F");
1968   m_pair_tree->Branch("pca_y", &m_pair.pca_y, "pca_y/F");
1969   m_pair_tree->Branch("pca_z", &m_pair.pca_z, "pca_z/F");
1970   m_pair_tree->Branch("pca1_x", &m_pair.pca1_x, "pca1_x/F");
1971   m_pair_tree->Branch("pca1_y", &m_pair.pca1_y, "pca1_y/F");
1972   m_pair_tree->Branch("pca1_z", &m_pair.pca1_z, "pca1_z/F");
1973   m_pair_tree->Branch("pca2_x", &m_pair.pca2_x, "pca2_x/F");
1974   m_pair_tree->Branch("pca2_y", &m_pair.pca2_y, "pca2_y/F");
1975   m_pair_tree->Branch("pca2_z", &m_pair.pca2_z, "pca2_z/F");
1976   m_pair_tree->Branch("v0_px", &m_pair.v0_px, "v0_px/F");
1977   m_pair_tree->Branch("v0_py", &m_pair.v0_py, "v0_py/F");
1978   m_pair_tree->Branch("v0_pz", &m_pair.v0_pz, "v0_pz/F");
1979   m_pair_tree->Branch("v0_pt", &m_pair.v0_pt, "v0_pt/F");
1980   m_pair_tree->Branch("mass_Kshort", &m_pair.mass_Kshort, "mass_Kshort/F");
1981   m_pair_tree->Branch("mass_Lambda", &m_pair.mass_Lambda, "mass_Lambda/F");
1982   m_pair_tree->Branch("mass_AntiLambda", &m_pair.mass_AntiLambda, "mass_AntiLambda/F");
1983   m_pair_tree->Branch("true_decay_x", &m_pair.true_decay_x, "true_decay_x/F");
1984   m_pair_tree->Branch("true_decay_y", &m_pair.true_decay_y, "true_decay_y/F");
1985   m_pair_tree->Branch("true_decay_z", &m_pair.true_decay_z, "true_decay_z/F");
1986   m_pair_tree->Branch("pca_to_true_3d", &m_pair.pca_to_true_3d, "pca_to_true_3d/F");
1987   m_pair_tree->Branch("pca_to_true_xy", &m_pair.pca_to_true_xy, "pca_to_true_xy/F");
1988   m_pair_tree->Branch("pca_to_true_z", &m_pair.pca_to_true_z, "pca_to_true_z/F");
1989   m_pair_tree->Branch("truth_alpha", &m_pair.truth_alpha, "truth_alpha/F");
1990   m_pair_tree->Branch("truth_qT", &m_pair.truth_qT, "truth_qT/F");
1991   m_pair_tree->Branch("delta_alpha", &m_pair.delta_alpha, "delta_alpha/F");
1992   m_pair_tree->Branch("delta_qT", &m_pair.delta_qT, "delta_qT/F");
1993   m_pair_tree->Branch("truth_px1", &m_pair.truth_px1, "truth_px1/F");
1994   m_pair_tree->Branch("truth_py1", &m_pair.truth_py1, "truth_py1/F");
1995   m_pair_tree->Branch("truth_pz1", &m_pair.truth_pz1, "truth_pz1/F");
1996   m_pair_tree->Branch("truth_px2", &m_pair.truth_px2, "truth_px2/F");
1997   m_pair_tree->Branch("truth_py2", &m_pair.truth_py2, "truth_py2/F");
1998   m_pair_tree->Branch("truth_pz2", &m_pair.truth_pz2, "truth_pz2/F");
1999   m_pair_tree->Branch("cos_mom1_truth", &m_pair.cos_mom1_truth, "cos_mom1_truth/F");
2000   m_pair_tree->Branch("cos_mom2_truth", &m_pair.cos_mom2_truth, "cos_mom2_truth/F");
2001   m_pair_tree->Branch("pca_theta1", &m_pair.pca_theta1, "pca_theta1/F");
2002   m_pair_tree->Branch("pca_theta2", &m_pair.pca_theta2, "pca_theta2/F");
2003   m_pair_tree->Branch("kalman_chi2_1", &m_pair.kalman_chi2_1, "kalman_chi2_1/F");
2004   m_pair_tree->Branch("kalman_chi2_2", &m_pair.kalman_chi2_2, "kalman_chi2_2/F");
2005   m_pair_tree->Branch("kalman_chi2_ndf1", &m_pair.kalman_chi2_ndf1, "kalman_chi2_ndf1/F");
2006   m_pair_tree->Branch("kalman_chi2_ndf2", &m_pair.kalman_chi2_ndf2, "kalman_chi2_ndf2/F");
2007   m_pair_tree->Branch("quality1", &m_pair.quality1, "quality1/F");
2008   m_pair_tree->Branch("quality2", &m_pair.quality2, "quality2/F");
2009   m_pair_tree->Branch("track_id1", &m_pair.track_id1, "track_id1/I");
2010   m_pair_tree->Branch("track_id2", &m_pair.track_id2, "track_id2/I");
2011   m_pair_tree->Branch("pid1", &m_pair.pid1, "pid1/I");
2012   m_pair_tree->Branch("pid2", &m_pair.pid2, "pid2/I");
2013   m_pair_tree->Branch("parent_id1", &m_pair.parent_id1, "parent_id1/I");
2014   m_pair_tree->Branch("parent_id2", &m_pair.parent_id2, "parent_id2/I");
2015   m_pair_tree->Branch("parent_pid", &m_pair.parent_pid, "parent_pid/I");
2016   m_pair_tree->Branch("kalman_ndof1", &m_pair.kalman_ndof1, "kalman_ndof1/I");
2017   m_pair_tree->Branch("kalman_ndof2", &m_pair.kalman_ndof2, "kalman_ndof2/I");
2018   m_pair_tree->Branch("npoints1", &m_pair.npoints1, "npoints1/S");
2019   m_pair_tree->Branch("npoints2", &m_pair.npoints2, "npoints2/S");
2020
2021   m_track_tree->Branch("run", &m_track.run, "run/I");
2022   m_track_tree->Branch("evt", &m_track.evt, "evt/I");
2023   m_track_tree->Branch("track_id", &m_track.track_id, "track_id/I");
2024   m_track_tree->Branch("shower_id", &m_track.shower_id, "shower_id/I");
2025   m_track_tree->Branch("pid", &m_track.pid, "pid/I");
2026   m_track_tree->Branch("parent_id", &m_track.parent_id, "parent_id/I");
2027   m_track_tree->Branch("parent_pid", &m_track.parent_pid, "parent_pid/I");
2028   m_track_tree->Branch("charge", &m_track.charge, "charge/D");
2029   m_track_tree->Branch("side", &m_track.side, "side/I");
2030   m_track_tree->Branch("npoints", &m_track.npoints, "npoints/I");
2031   m_track_tree->Branch("ntpc_clusters", &m_track.ntpc_clusters, "ntpc_clusters/i");
2032   m_track_tree->Branch("has_helix", &m_track.has_helix, "has_helix/I");
2033   m_track_tree->Branch("has_kalman", &m_track.has_kalman, "has_kalman/I");
2034   m_track_tree->Branch("is_primary", &m_track.is_primary, "is_primary/I");
2035   m_track_tree->Branch("px", &m_track.px, "px/D");
2036   m_track_tree->Branch("py", &m_track.py, "py/D");
2037   m_track_tree->Branch("pz", &m_track.pz, "pz/D");
2038   m_track_tree->Branch("pt", &m_track.pt, "pt/D");
2039   m_track_tree->Branch("p", &m_track.p, "p/D");
2040   m_track_tree->Branch("eta", &m_track.eta, "eta/D");
2041   m_track_tree->Branch("dedx", &m_track.dedx, "dedx/D");
2042   m_track_tree->Branch("x", &m_track.x, "x/F");
2043   m_track_tree->Branch("y", &m_track.y, "y/F");
2044   m_track_tree->Branch("z", &m_track.z, "z/F");
2045   m_track_tree->Branch("first_x", &m_track.first_x, "first_x/F");
2046   m_track_tree->Branch("first_y", &m_track.first_y, "first_y/F");
2047   m_track_tree->Branch("first_z", &m_track.first_z, "first_z/F");
2048   m_track_tree->Branch("first_r", &m_track.first_r, "first_r/F");
2049   m_track_tree->Branch("last_x", &m_track.last_x, "last_x/F");
2050   m_track_tree->Branch("last_y", &m_track.last_y, "last_y/F");
2051   m_track_tree->Branch("last_z", &m_track.last_z, "last_z/F");
2052   m_track_tree->Branch("last_r", &m_track.last_r, "last_r/F");
2053   m_track_tree->Branch("dca_xy", &m_track.dca_xy, "dca_xy/F");
2054   m_track_tree->Branch("dca_z", &m_track.dca_z, "dca_z/F");
2055   m_track_tree->Branch("vertex_x", &m_track.vertex_x, "vertex_x/D");
2056   m_track_tree->Branch("vertex_y", &m_track.vertex_y, "vertex_y/D");
2057   m_track_tree->Branch("vertex_z", &m_track.vertex_z, "vertex_z/D");
2058   m_track_tree->Branch("vertex_from_upstream", &m_track.vertex_from_upstream,
2059                        "vertex_from_upstream/I");
2060   m_track_tree->Branch("vertex_z_rms", &m_track.vertex_z_rms, "vertex_z_rms/D");
2061   m_track_tree->Branch("vertex_ntracks", &m_track.vertex_ntracks, "vertex_ntracks/i");
2062   m_track_tree->Branch("pca_x", &m_track.pca_x, "pca_x/D");
2063   m_track_tree->Branch("pca_y", &m_track.pca_y, "pca_y/D");
2064   m_track_tree->Branch("pca_z", &m_track.pca_z, "pca_z/D");
2065   m_track_tree->Branch("rDCA_zero", &m_track.rDCA_zero, "rDCA_zero/D");
2066   m_track_tree->Branch("zDCA", &m_track.zDCA, "zDCA/D");
2067   m_track_tree->Branch("helix_cx", &m_track.helix_cx, "helix_cx/F");
2068   m_track_tree->Branch("helix_cy", &m_track.helix_cy, "helix_cy/F");
2069   m_track_tree->Branch("helix_radius", &m_track.helix_radius, "helix_radius/F");
2070   m_track_tree->Branch("helix_z0", &m_track.helix_z0, "helix_z0/F");
2071   m_track_tree->Branch("helix_pitch", &m_track.helix_pitch, "helix_pitch/F");
2072   m_track_tree->Branch("helix_theta_first", &m_track.helix_theta_first, "helix_theta_first/F");
2073   m_track_tree->Branch("helix_theta_last", &m_track.helix_theta_last, "helix_theta_last/F");
2074   m_track_tree->Branch("helix_direction", &m_track.helix_direction, "helix_direction/F");
2075   m_track_tree->Branch("helix_search_anchored", &m_track.helix_search_anchored,
2076                        "helix_search_anchored/I");
2077   m_track_tree->Branch("helix_anchor_point_index", &m_track.helix_anchor_point_index,
2078                        "helix_anchor_point_index/I");
2079   m_track_tree->Branch("helix_anchor_theta", &m_track.helix_anchor_theta,
2080                        "helix_anchor_theta/F");
2081   m_track_tree->Branch("helix_anchor_path_cm", &m_track.helix_anchor_path_cm,
2082                        "helix_anchor_path_cm/F");
2083   m_track_tree->Branch("helix_anchor_residual_cm", &m_track.helix_anchor_residual_cm,
2084                        "helix_anchor_residual_cm/F");
2085   m_track_tree->Branch("helix_search_theta_min", &m_track.helix_search_theta_min,
2086                        "helix_search_theta_min/F");
2087   m_track_tree->Branch("helix_search_theta_max", &m_track.helix_search_theta_max,
2088                        "helix_search_theta_max/F");
2089   m_track_tree->Branch("helix_search_upstream_cm", &m_track.helix_search_upstream_cm,
2090                        "helix_search_upstream_cm/F");
2091   m_track_tree->Branch("helix_search_downstream_cm", &m_track.helix_search_downstream_cm,
2092                        "helix_search_downstream_cm/F");
2093   m_track_tree->Branch("kalman_chi2", &m_track.kalman_chi2, "kalman_chi2/F");
2094   m_track_tree->Branch("kalman_ndof", &m_track.kalman_ndof, "kalman_ndof/I");
2095   m_track_tree->Branch("kalman_naccepted", &m_track.kalman_naccepted, "kalman_naccepted/i");
2096   m_track_tree->Branch("kalman_nrejected", &m_track.kalman_nrejected, "kalman_nrejected/i");
2097   m_track_tree->Branch("kalman_measurement_chi2", &m_track.kalman_measurement_chi2);
2098   m_track_tree->Branch("kalman_measurement_used", &m_track.kalman_measurement_used);
2099   if (m_kalman_config.collect_innovation_components)
2100   {
2101     m_track_tree->Branch("kalman_measurement_sigma_r", &m_track.kalman_measurement_sigma_r,
2102                          "kalman_measurement_sigma_r/F");
2103     m_track_tree->Branch("kalman_measurement_sigma_rphi", &m_track.kalman_measurement_sigma_rphi,
2104                          "kalman_measurement_sigma_rphi/F");
2105     m_track_tree->Branch("kalman_measurement_sigma_z", &m_track.kalman_measurement_sigma_z,
2106                          "kalman_measurement_sigma_z/F");
2107     m_track_tree->Branch("kalman_measurement_in_seed", &m_track.kalman_measurement_in_seed);
2108     m_track_tree->Branch("kalman_innovation_residual_r", &m_track.kalman_innovation_residual_r);
2109     m_track_tree->Branch("kalman_innovation_residual_rphi", &m_track.kalman_innovation_residual_rphi);
2110     m_track_tree->Branch("kalman_innovation_residual_z", &m_track.kalman_innovation_residual_z);
2111     m_track_tree->Branch("kalman_prediction_sigma_r", &m_track.kalman_prediction_sigma_r);
2112     m_track_tree->Branch("kalman_prediction_sigma_rphi", &m_track.kalman_prediction_sigma_rphi);
2113     m_track_tree->Branch("kalman_prediction_sigma_z", &m_track.kalman_prediction_sigma_z);
2114     m_track_tree->Branch("kalman_innovation_sigma_r", &m_track.kalman_innovation_sigma_r);
2115     m_track_tree->Branch("kalman_innovation_sigma_rphi", &m_track.kalman_innovation_sigma_rphi);
2116     m_track_tree->Branch("kalman_innovation_sigma_z", &m_track.kalman_innovation_sigma_z);
2117     m_track_tree->Branch("kalman_innovation_rho_r_rphi", &m_track.kalman_innovation_rho_r_rphi);
2118     m_track_tree->Branch("kalman_innovation_rho_r_z", &m_track.kalman_innovation_rho_r_z);
2119     m_track_tree->Branch("kalman_innovation_rho_rphi_z", &m_track.kalman_innovation_rho_rphi_z);
2120     m_track_tree->Branch("kalman_innovation_whitened_0", &m_track.kalman_innovation_whitened_0);
2121     m_track_tree->Branch("kalman_innovation_whitened_1", &m_track.kalman_innovation_whitened_1);
2122     m_track_tree->Branch("kalman_innovation_whitened_2", &m_track.kalman_innovation_whitened_2);
2123   }
2124   m_track_tree->Branch("kalman_qop_t", &m_track.kalman_qop_t, "kalman_qop_t/F");
2125   m_track_tree->Branch("kalman_omega", &m_track.kalman_omega, "kalman_omega/F");
2126   m_track_tree->Branch("kalman_cx", &m_track.kalman_cx, "kalman_cx/F");
2127   m_track_tree->Branch("kalman_cy", &m_track.kalman_cy, "kalman_cy/F");
2128   m_track_tree->Branch("kalman_radius", &m_track.kalman_radius, "kalman_radius/F");
2129   m_track_tree->Branch("fit_chi2", &m_track.fit_chi2, "fit_chi2/F");
2130   m_track_tree->Branch("fit_ndf", &m_track.fit_ndf, "fit_ndf/I");
2131   m_track_tree->Branch("quality", &m_track.quality, "quality/F");
2132   m_track_tree->Branch("truth_px", &m_track.truth_px, "truth_px/F");
2133   m_track_tree->Branch("truth_py", &m_track.truth_py, "truth_py/F");
2134   m_track_tree->Branch("truth_pz", &m_track.truth_pz, "truth_pz/F");
2135   m_track_tree->Branch("cos_mom_truth", &m_track.cos_mom_truth, "cos_mom_truth/F");
2136   m_track_tree->Branch("cluster_index", &m_track.cluster_index);
2137   m_track_tree->Branch("cluster_side", &m_track.cluster_side);
2138   m_track_tree->Branch("layer", &m_track.layer);
2139   m_track_tree->Branch("cluster_z", &m_track.cluster_z);
2140   m_track_tree->Branch("cluster_r", &m_track.cluster_r);
2141   m_track_tree->Branch("cluster_phi", &m_track.cluster_phi);
2142   m_track_tree->Branch("residual_z", &m_track.residual_z);
2143   m_track_tree->Branch("residual_r", &m_track.residual_r);
2144   m_track_tree->Branch("residual_rphi", &m_track.residual_rphi);
2145
2146   if (!m_cluster_residual_tree)
2147   {
2148     return;
2149   }
2150
2151   m_cluster_residual_tree->Branch("run", &m_cluster_residual.run, "run/I");
2152   m_cluster_residual_tree->Branch("evt", &m_cluster_residual.evt, "evt/I");
2153   m_cluster_residual_tree->Branch("track_id", &m_cluster_residual.track_id, "track_id/I");
2154   m_cluster_residual_tree->Branch("charge", &m_cluster_residual.charge, "charge/I");
2155   m_cluster_residual_tree->Branch("side", &m_cluster_residual.side, "side/I");
2156   m_cluster_residual_tree->Branch("layer", &m_cluster_residual.layer, "layer/I");
2157   m_cluster_residual_tree->Branch("cluster_index", &m_cluster_residual.cluster_index, "cluster_index/I");
2158   m_cluster_residual_tree->Branch("ntp_cluster", &m_cluster_residual.ntp_cluster, "ntp_cluster/I");
2159   m_cluster_residual_tree->Branch("npoints", &m_cluster_residual.npoints, "npoints/I");
2160   m_cluster_residual_tree->Branch("has_helix", &m_cluster_residual.has_helix, "has_helix/I");
2161   m_cluster_residual_tree->Branch("has_kalman", &m_cluster_residual.has_kalman, "has_kalman/I");
2162   m_cluster_residual_tree->Branch("cluster_x", &m_cluster_residual.cluster_x, "cluster_x/F");
2163   m_cluster_residual_tree->Branch("cluster_y", &m_cluster_residual.cluster_y, "cluster_y/F");
2164   m_cluster_residual_tree->Branch("cluster_z", &m_cluster_residual.cluster_z, "cluster_z/F");
2165   m_cluster_residual_tree->Branch("cluster_r", &m_cluster_residual.cluster_r, "cluster_r/F");
2166   m_cluster_residual_tree->Branch("cluster_phi", &m_cluster_residual.cluster_phi, "cluster_phi/F");
2167   m_cluster_residual_tree->Branch("fit_x", &m_cluster_residual.fit_x, "fit_x/F");
2168   m_cluster_residual_tree->Branch("fit_y", &m_cluster_residual.fit_y, "fit_y/F");
2169   m_cluster_residual_tree->Branch("fit_z", &m_cluster_residual.fit_z, "fit_z/F");
2170   m_cluster_residual_tree->Branch("fit_r", &m_cluster_residual.fit_r, "fit_r/F");
2171   m_cluster_residual_tree->Branch("fit_phi", &m_cluster_residual.fit_phi, "fit_phi/F");
2172   m_cluster_residual_tree->Branch("residual_x", &m_cluster_residual.residual_x, "residual_x/F");
2173   m_cluster_residual_tree->Branch("residual_y", &m_cluster_residual.residual_y, "residual_y/F");
2174   m_cluster_residual_tree->Branch("residual_z", &m_cluster_residual.residual_z, "residual_z/F");
2175   m_cluster_residual_tree->Branch("residual_r", &m_cluster_residual.residual_r, "residual_r/F");
2176   m_cluster_residual_tree->Branch("residual_rphi", &m_cluster_residual.residual_rphi, "residual_rphi/F");
2177   m_cluster_residual_tree->Branch("fit_chi2", &m_cluster_residual.fit_chi2, "fit_chi2/F");
2178   m_cluster_residual_tree->Branch("fit_ndf", &m_cluster_residual.fit_ndf, "fit_ndf/I");
2179   m_cluster_residual_tree->Branch("fit_chi2_ndf", &m_cluster_residual.fit_chi2_ndf, "fit_chi2_ndf/F");
2180 }
2181
2182 int TpcV0CandidateTree::pdg_charge(const int pid)
2183 {
2184   const int apid = std::abs(pid);
2185   int charge = 0;
2186   switch (apid)
2187   {
2188   case 11:
2189   case 13:
2190     charge = -1;
2191     break;
2192   case 211:
2193   case 321:
2194   case 2212:
2195   case 3222:
2196     charge = 1;
2197     break;
2198   case 3112:
2199   case 3312:
2200   case 3334:
2201     charge = -1;
2202     break;
2203   default:
2204     charge = 0;
2205     break;
2206   }
2207   return (pid < 0) ? -charge : charge;
2208 }
2209
2210 float TpcV0CandidateTree::quiet_nan()
2211 {
2212   return std::numeric_limits<float>::quiet_NaN();
2213 }
2214
2215 bool TpcV0CandidateTree::parse_point_order(const std::string &mode, PointOrder &order)
2216 {
2217   return TpcTrackHelixFitter::parse_point_order(mode, order);
2218 }
2219
2220 bool TpcV0CandidateTree::finite(const Vec3 &value)
2221 {
2222   return TpcTrackHelixFitter::finite(value);
2223 }
2224
2225 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::add(const Vec3 &lhs, const Vec3 &rhs)
2226 {
2227   return TpcTrackHelixFitter::add(lhs, rhs);
2228 }
2229
2230 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::subtract(const Vec3 &lhs, const Vec3 &rhs)
2231 {
2232   return TpcTrackHelixFitter::subtract(lhs, rhs);
2233 }
2234
2235 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::scale(const Vec3 &value, const double factor)
2236 {
2237   return TpcTrackHelixFitter::scale(value, factor);
2238 }
2239
2240 double TpcV0CandidateTree::dot(const Vec3 &lhs, const Vec3 &rhs)
2241 {
2242   return TpcTrackHelixFitter::dot(lhs, rhs);
2243 }
2244
2245 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::cross(const Vec3 &lhs, const Vec3 &rhs)
2246 {
2247   return {
2248       lhs.y * rhs.z - lhs.z * rhs.y,
2249       lhs.z * rhs.x - lhs.x * rhs.z,
2250       lhs.x * rhs.y - lhs.y * rhs.x};
2251 }
2252
2253 double TpcV0CandidateTree::norm(const Vec3 &value)
2254 {
2255   return TpcTrackHelixFitter::norm(value);
2256 }
2257
2258 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::unit(const Vec3 &value)
2259 {
2260   return TpcTrackHelixFitter::unit(value);
2261 }
2262
2263 double TpcV0CandidateTree::pt(const Vec3 &value)
2264 {
2265   return TpcTrackHelixFitter::pt(value);
2266 }
2267
2268 double TpcV0CandidateTree::distance(const Vec3 &lhs, const Vec3 &rhs)
2269 {
2270   return TpcTrackHelixFitter::distance(lhs, rhs);
2271 }
2272
2273 double TpcV0CandidateTree::vector_cosine(const Vec3 &lhs, const Vec3 &rhs)
2274 {
2275   return TpcTrackHelixFitter::vector_cosine(lhs, rhs);
2276 }
2277
2278 bool TpcV0CandidateTree::fit_circle_least_squares(const std::vector<TruthPoint> &points,
2279                                                   const std::size_t nfit,
2280                                                   double &cx,
2281                                                   double &cy,
2282                                                   double &radius)
2283 {
2284   return TpcTrackHelixFitter::fit_circle_least_squares(points, nfit, cx, cy, radius);
2285 }
2286
2287 void TpcV0CandidateTree::order_track_points(std::vector<TruthPoint> &points, const PointOrder order)
2288 {
2289   TpcTrackHelixFitter::order_points(points, order);
2290 }
2291
2292 bool TpcV0CandidateTree::fit_helix(const std::vector<TruthPoint> &points,
2293                                    const int fit_first_points,
2294                                    const int charge,
2295                                    const double bfield_t,
2296                                    HelixFit &helix)
2297 {
2298   return TpcTrackHelixFitter::fit(points, fit_first_points, bfield_t, helix) &&
2299          TpcTrackHelixFitter::orient_to_charge(helix, charge);
2300 }
2301
2302 bool TpcV0CandidateTree::fit_kalman(const std::vector<TruthPoint> &points,
2303                                     const int charge,
2304                                     TpcKalmanResult &kalman) const
2305 {
2306   const auto start = std::chrono::steady_clock::now();
2307   const bool success = TpcTrackKalmanFitter::fit(
2308       points, charge, m_kalman_config, kalman, kPionMass);
2309   const double fit_seconds = std::chrono::duration<double>(
2310                                  std::chrono::steady_clock::now() - start)
2311                                  .count();
2312   m_timing_kalman_fit_seconds += fit_seconds;
2313   m_timing_rkn_seconds += kalman.rkn_seconds;
2314   m_timing_rkn_propagations += kalman.rkn_propagations;
2315   m_timing_rkn_accepted_steps += kalman.rkn_accepted_steps;
2316   m_timing_rkn_rejected_trials += kalman.rkn_rejected_trials;
2317   m_timing_rkn_failures += kalman.rkn_failures;
2318   ++m_timing_kalman_fits;
2319   if (m_print_timing &&
2320       (m_timing_kalman_fits <= 3 || m_timing_kalman_fits % 10 == 0 || fit_seconds > 1.0))
2321   {
2322     std::cout << "[V0TimingFit] fit=" << m_timing_kalman_fits
2323               << " points=" << points.size()
2324               << " success=" << success
2325               << " fit_s=" << fit_seconds
2326               << " rkn_s=" << kalman.rkn_seconds
2327               << " rkn_propagations=" << kalman.rkn_propagations
2328               << " rkn_steps=" << kalman.rkn_accepted_steps
2329               << " rkn_retries=" << kalman.rkn_rejected_trials
2330               << " rkn_failures=" << kalman.rkn_failures
2331               << std::endl;
2332   }
2333   return success;
2334 }
2335
2336 bool TpcV0CandidateTree::helix_from_state(const Vec3 &position, const Vec3 &momentum,
2337                                           const int charge, const double bfield_t,
2338                                           HelixFit &helix)
2339 {
2340   return TpcTrackHelixFitter::from_state(position, momentum, charge, bfield_t, helix);
2341 }
2342
2343 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::helix_point(const HelixFit &helix, const double theta)
2344 {
2345   return TpcTrackHelixFitter::point(helix, theta);
2346 }
2347
2348 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::helix_tangent(const HelixFit &helix, const double theta)
2349 {
2350   return TpcTrackHelixFitter::tangent(helix, theta);
2351 }
2352
2353 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::helix_momentum(const HelixFit &helix, const double theta)
2354 {
2355   return TpcTrackHelixFitter::momentum(helix, theta);
2356 }
2357
2358 std::pair<double, double> TpcV0CandidateTree::theta_search_range(const HelixFit &helix,
2359                                                                  const double theta_extension,
2360                                                                  const double downstream_margin)
2361 {
2362   return TpcTrackHelixFitter::theta_search_range(helix, theta_extension, downstream_margin);
2363 }
2364
2365 bool TpcV0CandidateTree::line_line_pca(const Vec3 &pos1, const Vec3 &dir1,
2366                                        const Vec3 &pos2, const Vec3 &dir2,
2367                                        LinePca &pca, const bool normalize_dirs)
2368 {
2369   return TpcTrackHelixFitter::line_line_pca(pos1, dir1, pos2, dir2, pca, normalize_dirs);
2370 }
2371
2372 TpcV0CandidateTree::HelixPca TpcV0CandidateTree::refine_helix_pair(
2373     const HelixFit &helix1, const HelixFit &helix2,
2374     double theta1, double theta2,
2375     const double min1, const double max1,
2376     const double min2, const double max2,
2377     double max_step)
2378 {
2379   return TpcTrackHelixFitter::refine_pair(helix1, helix2, theta1, theta2,
2380                                           min1, max1, min2, max2, max_step);
2381 }
2382
2383 std::vector<TpcV0CandidateTree::HelixPca> TpcV0CandidateTree::helix_helix_pca_candidates(
2384     const HelixFit &helix1, const HelixFit &helix2,
2385     const double theta_extension,
2386     const int coarse_steps,
2387     const double downstream_margin,
2388     const int max_candidates)
2389 {
2390   return TpcTrackHelixFitter::pca_candidates(helix1, helix2, theta_extension,
2391                                              coarse_steps, downstream_margin,
2392                                              max_candidates);
2393 }
2394
2395 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::kalman_point(const TpcKalmanResult &kalman,
2396                                                           const double s_cm,
2397                                                           const TpcKalmanConfig &config,
2398                                                           const Vec3 &reference_vertex)
2399 {
2400   if (!kalman.success || kalman.states_smoothed.empty())
2401   {
2402     const double nan = quiet_nan();
2403     return {nan, nan, nan};
2404   }
2405
2406   const auto state = TpcTrackKalmanFitter::propagate_state(
2407       TpcTrackKalmanFitter::propagation_state(kalman, reference_vertex), s_cm, config, kalman.mass_gev);
2408   return TpcTrackKalmanFitter::state_position(state);
2409 }
2410
2411 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::kalman_tangent(const TpcKalmanResult &kalman,
2412                                                             const double s_cm,
2413                                                             const TpcKalmanConfig &config,
2414                                                             const Vec3 &reference_vertex)
2415 {
2416   if (!kalman.success || kalman.states_smoothed.empty())
2417   {
2418     const double nan = quiet_nan();
2419     return {nan, nan, nan};
2420   }
2421
2422   const auto state = TpcTrackKalmanFitter::propagate_state(
2423       TpcTrackKalmanFitter::propagation_state(kalman, reference_vertex), s_cm, config, kalman.mass_gev);
2424   return TpcTrackKalmanFitter::state_tangent(state);
2425 }
2426
2427 TpcV0CandidateTree::Vec3 TpcV0CandidateTree::kalman_momentum(const TpcKalmanResult &kalman,
2428                                                              const double s_cm,
2429                                                              const TpcKalmanConfig &config,
2430                                                              const Vec3 &reference_vertex)
2431 {
2432   if (!kalman.success || kalman.states_smoothed.empty())
2433   {
2434     const double nan = quiet_nan();
2435     return {nan, nan, nan};
2436   }
2437
2438   const auto state = TpcTrackKalmanFitter::propagate_state(
2439       TpcTrackKalmanFitter::propagation_state(kalman, reference_vertex), s_cm, config, kalman.mass_gev);
2440   return TpcTrackKalmanFitter::state_momentum(state);
2441 }
2442
2443 TpcV0CandidateTree::KalmanPca TpcV0CandidateTree::refine_kalman_pair(
2444     const TpcKalmanResult &kalman1,
2445     const TpcKalmanResult &kalman2,
2446     const TpcKalmanConfig &config,
2447     const Vec3 &reference_vertex,
2448     double s1, double s2,
2449     const double min1, const double max1,
2450     const double min2, const double max2,
2451     double max_step,
2452     const int max_iterations)
2453 {
2454   KalmanPca best;
2455   best.s1 = s1;
2456   best.s2 = s2;
2457   best.pca1 = kalman_point(kalman1, s1, config, reference_vertex);
2458   best.pca2 = kalman_point(kalman2, s2, config, reference_vertex);
2459   double best_dca2 = square(distance(best.pca1, best.pca2));
2460
2461   for (int iter = 0; iter < std::max(1, max_iterations); ++iter)
2462   {
2463     const Vec3 pos1 = kalman_point(kalman1, s1, config, reference_vertex);
2464     const Vec3 pos2 = kalman_point(kalman2, s2, config, reference_vertex);
2465     const Vec3 tan1 = kalman_tangent(kalman1, s1, config, reference_vertex);
2466     const Vec3 tan2 = kalman_tangent(kalman2, s2, config, reference_vertex);
2467
2468     LinePca line_pca;
2469     if (!line_line_pca(pos1, tan1, pos2, tan2, line_pca, false))
2470     {
2471       break;
2472     }
2473
2474     const double step1 = std::clamp(line_pca.step1, -max_step, max_step);
2475     const double step2 = std::clamp(line_pca.step2, -max_step, max_step);
2476     if (std::abs(step1) < 1.0e-4 && std::abs(step2) < 1.0e-4)
2477     {
2478       break;
2479     }
2480
2481     const double candidate_s1 = std::clamp(s1 + step1, min1, max1);
2482     const double candidate_s2 = std::clamp(s2 + step2, min2, max2);
2483     const Vec3 candidate_pos1 = kalman_point(kalman1, candidate_s1, config, reference_vertex);
2484     const Vec3 candidate_pos2 = kalman_point(kalman2, candidate_s2, config, reference_vertex);
2485     const double candidate_dca2 = square(distance(candidate_pos1, candidate_pos2));
2486     if (candidate_dca2 < best_dca2)
2487     {
2488       s1 = candidate_s1;
2489       s2 = candidate_s2;
2490       best_dca2 = candidate_dca2;
2491     }
2492     else
2493     {
2494       max_step *= 0.5;
2495       if (max_step < 1.0e-3)
2496       {
2497         break;
2498       }
2499     }
2500   }
2501
2502   best.s1 = s1;
2503   best.s2 = s2;
2504   best.pca1 = kalman_point(kalman1, s1, config, reference_vertex);
2505   best.pca2 = kalman_point(kalman2, s2, config, reference_vertex);
2506   best.dca = distance(best.pca1, best.pca2);
2507   return best;
2508 }
2509
2510 std::vector<TpcV0CandidateTree::KalmanPca> TpcV0CandidateTree::kalman_pca_candidates(
2511     const TpcKalmanResult &kalman1,
2512     const TpcKalmanResult &kalman2,
2513     const TpcKalmanConfig &config,
2514     const Vec3 &reference_vertex,
2515     const double max_upstream_cm,
2516     const double downstream_margin_cm,
2517     const int coarse_steps,
2518     const int max_candidates)
2519 {
2520   std::vector<KalmanPca> candidates;
2521   if (!kalman1.success || !kalman2.success ||
2522       kalman1.states_smoothed.empty() || kalman2.states_smoothed.empty())
2523   {
2524     return candidates;
2525   }
2526
2527   const double min1 = -std::abs(max_upstream_cm);
2528   const double max1 = std::abs(downstream_margin_cm);
2529   const double min2 = -std::abs(max_upstream_cm);
2530   const double max2 = std::abs(downstream_margin_cm);
2531   const int nsteps = std::max(8, coarse_steps);
2532   const int ncandidates = std::max(1, max_candidates);
2533   const bool use_fast_field_pca =
2534       config.magnetic_field != nullptr && config.rkn_fast_field_pca;
2535   TpcKalmanConfig search_config = config;
2536   if (use_fast_field_pca)
2537   {
2538     // Use a cheap uniform-Bz trajectory only to locate promising path lengths.
2539     // The selected candidates are refined below with the full field map.
2540     search_config.magnetic_field = nullptr;
2541     search_config.rkn_step_tolerance = 0.0;
2542     search_config.rkn_max_step_cm = std::max(20.0, std::abs(config.rkn_max_step_cm));
2543   }
2544   const TpcKalmanConfig &coarse_config = use_fast_field_pca ? search_config : config;
2545
2546   std::vector<std::tuple<double, double, double>> seeds;
2547   const auto nsteps_size = static_cast<std::size_t>(nsteps);
2548   seeds.reserve(nsteps_size * nsteps_size);
2549   std::vector<double> s1_values(static_cast<std::size_t>(nsteps));
2550   std::vector<double> s2_values(static_cast<std::size_t>(nsteps));
2551   std::vector<Vec3> points1(static_cast<std::size_t>(nsteps));
2552   std::vector<Vec3> points2(static_cast<std::size_t>(nsteps));
2553   for (int i = 0; i < nsteps; ++i)
2554   {
2555     s1_values[static_cast<std::size_t>(i)] =
2556         min1 + (max1 - min1) * static_cast<double>(i) / static_cast<double>(nsteps - 1);
2557     s2_values[static_cast<std::size_t>(i)] =
2558         min2 + (max2 - min2) * static_cast<double>(i) / static_cast<double>(nsteps - 1);
2559     points1[static_cast<std::size_t>(i)] =
2560         kalman_point(kalman1, s1_values[static_cast<std::size_t>(i)], coarse_config, reference_vertex);
2561     points2[static_cast<std::size_t>(i)] =
2562         kalman_point(kalman2, s2_values[static_cast<std::size_t>(i)], coarse_config, reference_vertex);
2563   }
2564
2565   for (int i = 0; i < nsteps; ++i)
2566   {
2567     const double s1 = s1_values[static_cast<std::size_t>(i)];
2568     const Vec3 &pos1 = points1[static_cast<std::size_t>(i)];
2569     for (int j = 0; j < nsteps; ++j)
2570     {
2571       const double s2 = s2_values[static_cast<std::size_t>(j)];
2572       const Vec3 &pos2 = points2[static_cast<std::size_t>(j)];
2573       const double d2 = square(distance(pos1, pos2));
2574       if (std::isfinite(d2))
2575       {
2576         seeds.emplace_back(d2, s1, s2);
2577       }
2578     }
2579   }
2580
2581   std::sort(seeds.begin(), seeds.end(),
2582             [](const auto &lhs, const auto &rhs)
2583             { return std::get<0>(lhs) < std::get<0>(rhs); });
2584
2585   const double max_step = std::max(max1 - min1, max2 - min2) / static_cast<double>(nsteps);
2586   const int nrefine = std::min(ncandidates, static_cast<int>(seeds.size()));
2587   candidates.reserve(static_cast<std::size_t>(nrefine));
2588   for (int index = 0; index < nrefine; ++index)
2589   {
2590     double seed_s1 = std::get<1>(seeds[static_cast<std::size_t>(index)]);
2591     double seed_s2 = std::get<2>(seeds[static_cast<std::size_t>(index)]);
2592     if (use_fast_field_pca)
2593     {
2594       const auto surrogate = refine_kalman_pair(
2595           kalman1, kalman2, search_config, reference_vertex,
2596           seed_s1, seed_s2, min1, max1, min2, max2, max_step, 12);
2597       seed_s1 = surrogate.s1;
2598       seed_s2 = surrogate.s2;
2599     }
2600
2601     candidates.push_back(refine_kalman_pair(
2602         kalman1, kalman2, config, reference_vertex,
2603         seed_s1, seed_s2, min1, max1, min2, max2, max_step,
2604         use_fast_field_pca
2605             ? std::max(1, config.rkn_field_pca_refine_iterations)
2606             : 30));
2607   }
2608
2609   std::sort(candidates.begin(), candidates.end(),
2610             [](const KalmanPca &lhs, const KalmanPca &rhs)
2611             { return lhs.dca < rhs.dca; });
2612   return candidates;
2613 }
2614
2615 std::pair<double, double> TpcV0CandidateTree::track_dca_to_vertex(const Vec3 &pos,
2616                                                                   const Vec3 &mom,
2617                                                                   const Vec3 &vertex)
2618 {
2619   return TpcTrackHelixFitter::line_dca_to_vertex(pos, mom, vertex);
2620 }
2621
2622 std::pair<double, double> TpcV0CandidateTree::helix_dca_to_vertex(const HelixFit &helix,
2623                                                                   const Vec3 &vertex)
2624 {
2625   return TpcTrackHelixFitter::helix_dca_to_vertex(helix, vertex);
2626 }
2627
2628 std::pair<double, double> TpcV0CandidateTree::fitted_track_dca_to_vertex(
2629     const Tracklet &tracklet,
2630     const Vec3 &vertex) const
2631 {
2632   if (tracklet.has_kalman)
2633   {
2634     return TpcTrackKalmanFitter::dca_to_vertex(
2635         tracklet.kalman, vertex, &m_kalman_config);
2636   }
2637   if (tracklet.has_helix)
2638   {
2639     return helix_dca_to_vertex(tracklet.helix, vertex);
2640   }
2641   return track_dca_to_vertex(tracklet.position, tracklet.momentum, vertex);
2642 }
2643
2644 bool TpcV0CandidateTree::armenteros(const Vec3 &pplus, const Vec3 &pminus,
2645                                     double &alpha, double &qt)
2646 {
2647   const Vec3 v0p = add(pplus, pminus);
2648   const Vec3 direction = unit(v0p);
2649   if (!finite(direction))
2650   {
2651     return false;
2652   }
2653
2654   const double pl_plus = dot(pplus, direction);
2655   const double pl_minus = dot(pminus, direction);
2656   const double denom = pl_plus + pl_minus;
2657   if (std::abs(denom) < 1e-10)
2658   {
2659     return false;
2660   }
2661
2662   alpha = (pl_plus - pl_minus) / denom;
2663   qt = norm(subtract(pplus, scale(direction, pl_plus)));
2664   return true;
2665 }
2666
2667 double TpcV0CandidateTree::invariant_mass(const Vec3 &mom1, const double mass1,
2668                                           const Vec3 &mom2, const double mass2)
2669 {
2670   const double e1 = std::sqrt(dot(mom1, mom1) + square(mass1));
2671   const double e2 = std::sqrt(dot(mom2, mom2) + square(mass2));
2672   const Vec3 total_mom = add(mom1, mom2);
2673   const double mass2_total = square(e1 + e2) - dot(total_mom, total_mom);
2674   return (mass2_total > 0.0) ? std::sqrt(mass2_total) : 0.0;
2675 }
2676
2677 bool TpcV0CandidateTree::passes_preselection(const Tracklet &track1, const Tracklet &track2,
2678                                              const Vec3 &primary_vertex) const
2679 {
2680   if (m_pre_track_npoints_min > 0 &&
2681       (track1.npoints < m_pre_track_npoints_min || track2.npoints < m_pre_track_npoints_min))
2682   {
2683     return false;
2684   }
2685
2686   if (m_pre_track_quality_max >= 0.0 &&
2687       (!std::isfinite(track1.fit_chi2_ndf) || !std::isfinite(track2.fit_chi2_ndf) ||
2688        track1.fit_chi2_ndf >= m_pre_track_quality_max ||
2689        track2.fit_chi2_ndf >= m_pre_track_quality_max))
2690   {
2691     return false;
2692   }
2693
2694   if (m_pre_track_pt_min > 0.0 &&
2695       (pt(track1.momentum) < m_pre_track_pt_min || pt(track2.momentum) < m_pre_track_pt_min))
2696   {
2697     return false;
2698   }
2699
2700   auto dca1 = track1.vertex_dca;
2701   if (!track1.has_vertex_dca)
2702   {
2703     dca1 = fitted_track_dca_to_vertex(track1, primary_vertex);
2704   }
2705   auto dca2 = track2.vertex_dca;
2706   if (!track2.has_vertex_dca)
2707   {
2708     dca2 = fitted_track_dca_to_vertex(track2, primary_vertex);
2709   }
2710   if (!std::isfinite(dca1.first) || !std::isfinite(dca1.second) ||
2711       !std::isfinite(dca2.first) || !std::isfinite(dca2.second))
2712   {
2713     return false;
2714   }
2715
2716   if (m_pre_track_dca_xy_min >= 0.0 &&
2717       (dca1.first < m_pre_track_dca_xy_min || dca2.first < m_pre_track_dca_xy_min))
2718   {
2719     return false;
2720   }
2721   if (m_pre_track_dca_z_min >= 0.0 &&
2722       (dca1.second < m_pre_track_dca_z_min || dca2.second < m_pre_track_dca_z_min))
2723   {
2724     return false;
2725   }
2726   if (m_pre_track_dca_xy_max >= 0.0 &&
2727       (dca1.first > m_pre_track_dca_xy_max || dca2.first > m_pre_track_dca_xy_max))
2728   {
2729     return false;
2730   }
2731   if (m_pre_track_dca_z_max >= 0.0 &&
2732       (dca1.second > m_pre_track_dca_z_max || dca2.second > m_pre_track_dca_z_max))
2733   {
2734     return false;
2735   }
2736
2737   if (m_pre_pair_dca_max < 0.0 && m_pre_lproj_min < 0.0 && m_pre_cos_theta_min < -1.0)
2738   {
2739     return true;
2740   }
2741
2742   LinePca rough_pca;
2743   if (!line_line_pca(track1.position, track1.momentum, track2.position, track2.momentum, rough_pca, true))
2744   {
2745     return false;
2746   }
2747
2748   if (m_pre_pair_dca_max >= 0.0 && rough_pca.dca > m_pre_pair_dca_max)
2749   {
2750     return false;
2751   }
2752
2753   const Vec3 rough_vertex = scale(add(rough_pca.pca1, rough_pca.pca2), 0.5);
2754   const Vec3 flight = subtract(rough_vertex, primary_vertex);
2755   const Vec3 total_mom = add(track1.momentum, track2.momentum);
2756   const double lproj = norm(flight);
2757   const double cos_theta = vector_cosine(flight, total_mom);
2758
2759   if (m_pre_lproj_min >= 0.0 && (!std::isfinite(lproj) || lproj < m_pre_lproj_min))
2760   {
2761     return false;
2762   }
2763   if (m_pre_cos_theta_min >= -1.0 &&
2764       (!std::isfinite(cos_theta) || cos_theta < m_pre_cos_theta_min))
2765   {
2766     return false;
2767   }
2768
2769   return true;
2770 }
2771
2772 bool TpcV0CandidateTree::passes_pair_selection(const Vec3 &pca1, const Vec3 &pca2,
2773                                                const Vec3 &pair_vertex,
2774                                                const Vec3 &primary_vertex,
2775                                                const double pair_dca,
2776                                                const double cos_theta,
2777                                                const double alpha) const
2778 {
2779   if (m_pair_pca_z_max >= 0.0 &&
2780       (!std::isfinite(pair_vertex.z) || std::abs(pair_vertex.z) >= m_pair_pca_z_max))
2781   {
2782     return false;
2783   }
2784
2785   const double pca_dz = std::abs(pca1.z - pca2.z);
2786   if (m_pair_pca_dz_max >= 0.0 &&
2787       (!std::isfinite(pca_dz) || pca_dz >= m_pair_pca_dz_max))
2788   {
2789     return false;
2790   }
2791
2792   const double dx = pair_vertex.x - primary_vertex.x;
2793   const double dy = pair_vertex.y - primary_vertex.y;
2794   const double decay_radius = std::hypot(dx, dy);
2795   if (m_pair_decay_radius_min >= 0.0 &&
2796       (!std::isfinite(decay_radius) || decay_radius <= m_pair_decay_radius_min))
2797   {
2798     return false;
2799   }
2800
2801   if (m_pair_alpha_abs_max >= 0.0 &&
2802       (!std::isfinite(alpha) || std::abs(alpha) >= m_pair_alpha_abs_max))
2803   {
2804     return false;
2805   }
2806
2807   if (m_pair_dca_max >= 0.0 &&
2808       (!std::isfinite(pair_dca) || pair_dca >= m_pair_dca_max))
2809   {
2810     return false;
2811   }
2812
2813   if (m_pair_dira_min >= -1.0 &&
2814       (!std::isfinite(cos_theta) || cos_theta <= m_pair_dira_min))
2815   {
2816     return false;
2817   }
2818
2819   return true;
2820 }