packages/tpctrackreco/Tpc_AssembledTrackReco.cc

0001 #include "Tpc_AssembledTrackReco.h"
0002
0003 #include "IdealPadMap.h"
0004 #include "Tpc_AssembledTrack.h"
0005 #include "Tpc_AssembledTrackContainer.h"
0006 #include "Tpc_AssembledTrackContainerv1.h"
0007 #include "Tpc_AssembledTrackv1.h"
0008 #include "Tpc_FittingTools.h"
0009
0010 #include "Tpc_ModuleTrack.h"
0011 #include "Tpc_ModuleTrackContainer.h"
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/PHIODataNode.h>
0017 #include <phool/PHNodeIterator.h>
0018 #include <phool/PHObject.h>
0019 #include <phool/getClass.h>
0020
0021 #include <trackbase/TpcDefs.h>
0022 #include <trackbase/TrkrDefs.h>
0023 #include <trackbase/TrkrHit.h>
0024 #include <trackbase/TrkrHitSet.h>
0025 #include <trackbase/TrkrHitSetContainer.h>
0026
0027 #include <TFile.h>
0028 #include <TH1D.h>
0029 #include <TH2D.h>
0030 #include <TTree.h>
0031
0032 #include <algorithm>
0033 #include <cmath>
0034 #include <format>
0035 #include <functional>
0036 #include <iostream>
0037 #include <limits>
0038 #include <thread>
0039 #include <utility>
0040 #include <vector>
0041
0042 namespace
0043 {
0044   // NOLINTBEGIN(misc-non-private-member-variables-in-classes)
0045   struct PieceStartSort
0046   {
0047     const std::vector<Tpc_AssembledTrackReco::Piece>* pieces;
0048     explicit PieceStartSort(const std::vector<Tpc_AssembledTrackReco::Piece>* p)
0049       : pieces(p)
0050     {
0051     }
0052     bool operator()(unsigned int a, unsigned int b) const
0053     {
0054       const Tpc_AssembledTrackReco::Piece& pa = (*pieces)[a];
0055       const Tpc_AssembledTrackReco::Piece& pb = (*pieces)[b];
0056       if (pa.first_layer != pb.first_layer)
0057       {
0058         return pa.first_layer < pb.first_layer;
0059       }
0060       if (pa.last_layer != pb.last_layer)
0061       {
0062         return pa.last_layer < pb.last_layer;
0063       }
0064       if (pa.sector != pb.sector)
0065       {
0066         return pa.sector < pb.sector;
0067       }
0068       return pa.source_track_id < pb.source_track_id;
0069     }
0070   };
0071
0072   struct CandidateStartSort
0073   {
0074     const std::vector<Tpc_AssembledTrackReco::Candidate>* candidates;
0075     explicit CandidateStartSort(const std::vector<Tpc_AssembledTrackReco::Candidate>* c)
0076       : candidates(c)
0077     {
0078     }
0079     bool operator()(unsigned int a, unsigned int b) const
0080     {
0081       const Tpc_AssembledTrackReco::Candidate& ca = (*candidates)[a];
0082       const Tpc_AssembledTrackReco::Candidate& cb = (*candidates)[b];
0083       if (ca.first_layer != cb.first_layer)
0084       {
0085         return ca.first_layer < cb.first_layer;
0086       }
0087       if (ca.last_layer != cb.last_layer)
0088       {
0089         return ca.last_layer < cb.last_layer;
0090       }
0091       if (ca.first_sector != cb.first_sector)
0092       {
0093         return ca.first_sector < cb.first_sector;
0094       }
0095       return ca.nsegments < cb.nsegments;
0096     }
0097   };
0098
0099   struct RadiusSort
0100   {
0101     const std::vector<double>* radius;
0102     explicit RadiusSort(const std::vector<double>* r)
0103       : radius(r)
0104     {
0105     }
0106     bool operator()(unsigned int a, unsigned int b) const { return (*radius)[a] < (*radius)[b]; }
0107   };
0108   // NOLINTEND(misc-non-private-member-variables-in-classes)
0109
0110   double unwrap_phi_to_reference(double phi, const double ref)
0111   {
0112     while (phi - ref > M_PI)
0113     {
0114       phi -= 2.0 * M_PI;
0115     }
0116     while (phi - ref < -M_PI)
0117     {
0118       phi += 2.0 * M_PI;
0119     }
0120     return phi;
0121   }
0122
0123   double wrap_to_pi(double phi)
0124   {
0125     while (phi > M_PI)
0126     {
0127       phi -= 2.0 * M_PI;
0128     }
0129     while (phi <= -M_PI)
0130     {
0131       phi += 2.0 * M_PI;
0132     }
0133     return phi;
0134   }
0135
0136   int wrapped_sector_delta(const unsigned int sector_a, const unsigned int sector_b)
0137   {
0138     int d = static_cast<int>(sector_b) - static_cast<int>(sector_a);
0139     while (d > 6)
0140     {
0141       d -= 12;
0142     }
0143     while (d < -6)
0144     {
0145       d += 12;
0146     }
0147     return d;
0148   }
0149
0150   double sagitta_model_derivative(double xrot, double x0, double invR)
0151   {
0152     const double dx = xrot - x0;
0153     const double dx2 = dx * dx;
0154     const double invR2 = invR * invR;
0155     const double invR3 = invR2 * invR;
0156     const double invR5 = invR3 * invR2;
0157     return -invR * dx - 0.5 * invR3 * dx2 * dx - 0.375 * invR5 * dx2 * dx2 * dx;
0158   }
0159
0160   double predict_sagitta_phi(double radius, double S, double x0, double invR, double theta, double bline)
0161   {
0162     const double c = std::cos(theta);
0163     const double s = std::sin(theta);
0164     double yy = std::tan(theta) * radius;
0165
0166     for (unsigned int iter = 0; iter < 25; ++iter)
0167     {
0168       const double xrot = c * radius + s * yy;
0169       const double yrot = -s * radius + c * yy;
0170       const double f = Tpc_FittingTools::sagittaModel(xrot, S, x0, invR);
0171       const double g = yrot - f;
0172       const double df = sagitta_model_derivative(xrot, x0, invR);
0173       const double dg = c - df * s;
0174       if (std::fabs(dg) < 1.0e-12)
0175       {
0176         break;
0177       }
0178       const double step = g / dg;
0179       yy -= step;
0180       if (std::fabs(step) < 1.0e-10)
0181       {
0182         break;
0183       }
0184     }
0185
0186     return bline + yy;
0187   }
0188
0189   bool fit_points(const std::vector<double>& radius,
0190                   const std::vector<double>& phi,
0191                   const std::vector<double>& tbin,
0192                   const std::vector<double>& weight,
0193                   bool use_sagitta,
0194                   double& phi_slope,
0195                   double& phi_intercept,
0196                   double& phi_S,
0197                   double& phi_x0,
0198                   double& phi_invR,
0199                   double& phi_theta,
0200                   double& phi_bline,
0201                   bool& phi_sagitta_ok,
0202                   double& tbin_slope,
0203                   double& tbin_intercept,
0204                   double& chi2_phi,
0205                   double& chi2_tbin,
0206                   int& ndof_phi,
0207                   int& ndof_tbin)
0208   {
0209     if (radius.size() < 2 || radius.size() != phi.size() || radius.size() != tbin.size() || radius.size() != weight.size())
0210     {
0211       return false;
0212     }
0213
0214     std::vector<Tpc_FittingTools::FitPoint> phi_points;
0215     std::vector<Tpc_FittingTools::FitPoint> tbin_points;
0216     phi_points.reserve(radius.size());
0217     tbin_points.reserve(radius.size());
0218     for (unsigned int i = 0; i < radius.size(); ++i)
0219     {
0220       phi_points.emplace_back(radius[i], phi[i], weight[i]);
0221       tbin_points.emplace_back(radius[i], tbin[i], weight[i]);
0222     }
0223
0224     const Tpc_FittingTools::LineFit phi_line = Tpc_FittingTools::fitLine(phi_points);
0225     if (!phi_line.ok)
0226     {
0227       return false;
0228     }
0229
0230     phi_slope = phi_line.slope;
0231     phi_intercept = phi_line.intercept;
0232     chi2_phi = phi_line.chi2;
0233     ndof_phi = phi_line.ndof;
0234     phi_S = 0.0;
0235     phi_x0 = 0.0;
0236     phi_invR = 0.0;
0237     phi_theta = std::atan(phi_slope);
0238     phi_bline = phi_intercept;
0239     phi_sagitta_ok = false;
0240
0241     if (use_sagitta && radius.size() >= 3)
0242     {
0243       const Tpc_FittingTools::SagittaFit phi_sagitta = Tpc_FittingTools::fitSagitta(phi_points);
0244       if (phi_sagitta.ok)
0245       {
0246         phi_S = phi_sagitta.S;
0247         phi_x0 = phi_sagitta.x0;
0248         phi_invR = phi_sagitta.invR;
0249         phi_theta = phi_sagitta.theta;
0250         phi_bline = phi_sagitta.b;
0251         chi2_phi = phi_sagitta.chi2;
0252         ndof_phi = phi_sagitta.ndof;
0253         phi_sagitta_ok = true;
0254       }
0255     }
0256
0257     const Tpc_FittingTools::LineFit tbin_line = Tpc_FittingTools::fitLine(tbin_points);
0258     if (!tbin_line.ok)
0259     {
0260       return false;
0261     }
0262
0263     tbin_slope = tbin_line.slope;
0264     tbin_intercept = tbin_line.intercept;
0265     chi2_tbin = tbin_line.chi2;
0266     ndof_tbin = tbin_line.ndof;
0267     return true;
0268   }
0269 }  // namespace
0270
0271 Tpc_AssembledTrackReco::Piece::Piece()
0272   : source_index(0)
0273   , source_track_id(0)
0274   , event(0)
0275   , region(0)
0276   , sector(0)
0277   , side(0)
0278   , first_layer(0)
0279   , last_layer(0)
0280   , nblobs(0)
0281   , nrawhits(0)
0282   , phi_slope(0.0)
0283   , phi_intercept(0.0)
0284   , phi_S(0.0)
0285   , phi_x0(0.0)
0286   , phi_invR(0.0)
0287   , phi_theta(0.0)
0288   , phi_bline(0.0)
0289   , phi_sagitta_ok(false)
0290   , tbin_slope(0.0)
0291   , tbin_intercept(0.0)
0292 {
0293 }
0294
0295 Tpc_AssembledTrackReco::Candidate::Candidate()
0296   : event(0)
0297   , side(0)
0298   , first_layer(0)
0299   , last_layer(0)
0300   , first_sector(0)
0301   , last_sector(0)
0302   , first_region(0)
0303   , last_region(0)
0304   , nsegments(0)
0305   , nblobs(0)
0306   , nrawhits(0)
0307   , phi_slope(0.0)
0308   , phi_intercept(0.0)
0309   , phi_S(0.0)
0310   , phi_x0(0.0)
0311   , phi_invR(0.0)
0312   , phi_theta(0.0)
0313   , phi_bline(0.0)
0314   , phi_sagitta_ok(false)
0315   , tbin_slope_r(0.0)
0316   , tbin_intercept_r(0.0)
0317   , chi2_phi(0.0)
0318   , chi2_tbin(0.0)
0319   , ndof_phi(0)
0320   , ndof_tbin(0)
0321 {
0322 }
0323
0324 Tpc_AssembledTrackReco::Tpc_AssembledTrackReco(const std::string& name, const std::string& filename)
0325   : SubsysReco(name)
0326   , m_outputFileName(filename)
0327   , m_debugOutputFileName("Tpc_AssembledTrackRecoDebug.root")
0328   , m_inputNodeName("TPC_MODULETRACKS")
0329   , m_outputNodeName("TPC_ASSEMBLEDTRACKS")
0330   , m_outputFile(nullptr)
0331   , m_debugOutputFile(nullptr)
0332   , m_tree(nullptr)
0333   , m_tpcModuleTrackContainer(nullptr)
0334   , m_assembledTrackContainer(nullptr)
0335   , m_hits(nullptr)
0336   , m_event(0)
0337   , m_idealPadMap(nullptr)
0338   , m_connectMaxLayerGap(16)
0339   , m_connect_dphi(0.03)
0340   , m_connect_dtbin(8.0)
0341   , m_connect_dphi_slope(0.01)
0342   , m_connect_dtbin_slope(2.0)
0343   , m_useSagittaPhiFit(true)
0344   , m_seedSigmaX(5.0)
0345   , m_seedSigmaY(5.0)
0346   , m_seedSigmaZ(10.0)
0347   , m_seedSigmaPx(1.0)
0348   , m_seedSigmaPy(1.0)
0349   , m_seedSigmaPz(1.0)
0350   , m_h_dphi(nullptr)
0351   , m_h_dtbin(nullptr)
0352   , m_h_dmphi(nullptr)
0353   , m_h_dmtbin(nullptr)
0354   , m_h_score(nullptr)
0355   , m_h_dphi_vs_dtbin(nullptr)
0356   , m_h_dmphi_vs_dmtbin(nullptr)
0357   , m_h_dphi_vs_dmphi(nullptr)
0358   , m_h_tbin_slope_vs_first_tbin(nullptr)
0359   , m_h_tbin_slope_vs_last_tbin(nullptr)
0360   , m_h_track_tbin_slope_vs_tbin_span_3modules(nullptr)
0361   , m_h_track_tbin_slope_vs_first_tbin_3modules(nullptr)
0362   , m_h_track_tbin_slope_vs_last_tbin_3modules(nullptr)
0363   , m_h_layer_gap(nullptr)
0364   , m_h_nsegments(nullptr)
0365   , m_h_matched_sector_delta(nullptr)
0366 {
0367 }
0368
0369 Tpc_AssembledTrackReco::~Tpc_AssembledTrackReco()
0370 {
0371   delete m_idealPadMap;
0372   m_idealPadMap = nullptr;
0373
0374   if (m_debugOutputFile)
0375   {
0376     delete m_debugOutputFile;
0377     m_debugOutputFile = nullptr;
0378   }
0379
0380   if (m_outputFile)
0381   {
0382     delete m_outputFile;
0383     m_outputFile = nullptr;
0384   }
0385 }
0386
0387 int Tpc_AssembledTrackReco::Init(PHCompositeNode* /*unused*/)
0388 {
0389   if (Verbosity() <= 0)
0390   {
0391     return Fun4AllReturnCodes::EVENT_OK;
0392   }
0393
0394   m_outputFile = new TFile(m_outputFileName.c_str(), "RECREATE");
0395   if (!m_outputFile || m_outputFile->IsZombie())
0396   {
0397     std::cerr << Name() << "::Init - cannot create " << m_outputFileName << std::endl;
0398     return Fun4AllReturnCodes::ABORTRUN;
0399   }
0400
0401   m_debugOutputFile = new TFile(m_debugOutputFileName.c_str(), "RECREATE");
0402   if (!m_debugOutputFile || m_debugOutputFile->IsZombie())
0403   {
0404     std::cerr << Name() << "::Init - cannot create debug file " << m_debugOutputFileName << std::endl;
0405     return Fun4AllReturnCodes::ABORTRUN;
0406   }
0407   create_debug_histograms();
0408
0409   m_tree = new TTree("Tpc_AssembledTracks", "Assembled tracks connected from Tpc_ModuleTrackReco");
0410   m_tree->Branch("event", &m_tree_event, "event/I");
0411   m_tree->Branch("track_id", &m_tree_track_id);
0412   m_tree->Branch("side", &m_tree_side);
0413   m_tree->Branch("nsegments", &m_tree_nsegments);
0414   m_tree->Branch("nblobs", &m_tree_nblobs);
0415   m_tree->Branch("nrawhits", &m_tree_nrawhits);
0416   m_tree->Branch("first_layer", &m_tree_first_layer);
0417   m_tree->Branch("last_layer", &m_tree_last_layer);
0418   m_tree->Branch("first_sector", &m_tree_first_sector);
0419   m_tree->Branch("last_sector", &m_tree_last_sector);
0420   m_tree->Branch("first_region", &m_tree_first_region);
0421   m_tree->Branch("last_region", &m_tree_last_region);
0422
0423   m_tree->Branch("source_assembled_track_id", &m_tree_source_assembled_track_id);
0424   m_tree->Branch("source_inmodule_track_id", &m_tree_source_inmodule_track_id);
0425   m_tree->Branch("source_region", &m_tree_source_region);
0426   m_tree->Branch("source_sector", &m_tree_source_sector);
0427   m_tree->Branch("source_side", &m_tree_source_side);
0428
0429   m_tree->Branch("hit_assembled_track_id", &m_tree_hit_assembled_track_id);
0430   m_tree->Branch("hit_hitsetkey", &m_tree_hit_hitsetkey);
0431   m_tree->Branch("hit_hitkey", &m_tree_hit_hitkey);
0432
0433   return Fun4AllReturnCodes::EVENT_OK;
0434 }
0435
0436 int Tpc_AssembledTrackReco::InitRun(PHCompositeNode* topNode)
0437 {
0438   if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0439   {
0440     return Fun4AllReturnCodes::ABORTRUN;
0441   }
0442   if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0443   {
0444     return Fun4AllReturnCodes::ABORTRUN;
0445   }
0446
0447   delete m_idealPadMap;
0448   m_idealPadMap = new IdealPadMap();
0449   if (m_idealPadMap->load_from_cdb(Verbosity()) != 0 || !m_idealPadMap->is_loaded())
0450   {
0451     std::cerr << Name() << "::InitRun - cannot load IdealPadMap from CDB" << std::endl;
0452     return Fun4AllReturnCodes::ABORTRUN;
0453   }
0454
0455   m_event = 0;
0456   return Fun4AllReturnCodes::EVENT_OK;
0457 }
0458
0459 int Tpc_AssembledTrackReco::End(PHCompositeNode* /*unused*/)
0460 {
0461   if (m_outputFile)
0462   {
0463     m_outputFile->cd();
0464     if (m_tree)
0465     {
0466       m_tree->Write();
0467     }
0468     m_outputFile->Close();
0469     delete m_outputFile;
0470     m_outputFile = nullptr;
0471   }
0472
0473   write_debug_histograms();
0474   if (m_debugOutputFile)
0475   {
0476     m_debugOutputFile->Close();
0477     delete m_debugOutputFile;
0478     m_debugOutputFile = nullptr;
0479   }
0480
0481   return Fun4AllReturnCodes::EVENT_OK;
0482 }
0483
0484 int Tpc_AssembledTrackReco::getNodes(PHCompositeNode* topNode)
0485 {
0486   m_tpcModuleTrackContainer = findNode::getClass<Tpc_ModuleTrackContainer>(topNode, m_inputNodeName);
0487   if (!m_tpcModuleTrackContainer)
0488   {
0489     std::cerr << Name() << "::getNodes - missing " << m_inputNodeName << std::endl;
0490     return Fun4AllReturnCodes::ABORTRUN;
0491   }
0492
0493   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0494   if (!m_hits)
0495   {
0496     std::cerr << Name() << "::getNodes - missing TRKR_HITSET" << std::endl;
0497     return Fun4AllReturnCodes::ABORTRUN;
0498   }
0499
0500   return Fun4AllReturnCodes::EVENT_OK;
0501 }
0502
0503 int Tpc_AssembledTrackReco::createNodes(PHCompositeNode* topNode)
0504 {
0505   PHNodeIterator iter(topNode);
0506   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0507
0508   if (!dstNode)
0509   {
0510     dstNode = new PHCompositeNode("DST");
0511     topNode->addNode(dstNode);
0512   }
0513
0514   m_assembledTrackContainer = findNode::getClass<Tpc_AssembledTrackContainer>(topNode, m_outputNodeName);
0515   if (!m_assembledTrackContainer)
0516   {
0517     m_assembledTrackContainer = new Tpc_AssembledTrackContainerv1();
0518     PHIODataNode<PHObject>* node = new PHIODataNode<PHObject>(m_assembledTrackContainer, m_outputNodeName, "PHObject");
0519     dstNode->addNode(node);
0520     std::cout << Name() << "::createNodes - created " << m_outputNodeName << " node" << std::endl;
0521   }
0522
0523   return Fun4AllReturnCodes::EVENT_OK;
0524 }
0525
0526 void Tpc_AssembledTrackReco::create_debug_histograms()
0527 {
0528   if (!m_debugOutputFile)
0529   {
0530     return;
0531   }
0532   m_debugOutputFile->cd();
0533
0534   m_h_dphi = new TH1D("h_dphi", std::format("#Delta#phi at match point, cut={:.4g};|#Delta#phi| [rad];tested pairs", m_connect_dphi).c_str(), 200, 0.0, std::max(0.2, 5.0 * m_connect_dphi));
0535   m_h_dtbin = new TH1D("h_dtbin", std::format("#Deltatbin at match point, cut={:.4g};|#Deltatbin|;tested pairs", m_connect_dtbin).c_str(), 200, 0.0, std::max(50.0, 5.0 * m_connect_dtbin));
0536   m_h_dmphi = new TH1D("h_dmphi", std::format("#Delta(d#phi/dr), cut={:.4g};|#Delta(d#phi/dr)| [rad/cm];tested pairs", m_connect_dphi_slope).c_str(), 200, 0.0, std::max(0.08, 5.0 * m_connect_dphi_slope));
0537   m_h_dmtbin = new TH1D("h_dmtbin", std::format("#Delta(dtbin/dr), cut={:.4g};|#Delta(dtbin/dr)| [tbin/cm];tested pairs", m_connect_dtbin_slope).c_str(), 200, 0.0, std::max(20.0, 5.0 * m_connect_dtbin_slope));
0538   m_h_score = new TH1D("h_score", "accepted connection score;score;accepted connections", 200, 0.0, 20.0);
0539
0540   m_h_dphi_vs_dtbin = new TH2D("h_dphi_vs_dtbin", std::format("#Delta#phi vs #Deltatbin, cuts #Delta#phi<{:.4g} #Deltatbin<{:.4g};|#Delta#phi| [rad];|#Deltatbin|", m_connect_dphi, m_connect_dtbin).c_str(), 160, 0.0, std::max(0.2, 5.0 * m_connect_dphi), 160, 0.0, std::max(50.0, 5.0 * m_connect_dtbin));
0541   m_h_dmphi_vs_dmtbin = new TH2D("h_dmphi_vs_dmtbin", std::format("slope residuals, cuts #Delta(d#phi/dr)<{:.4g} #Delta(dtbin/dr)<{:.4g};|#Delta(d#phi/dr)| [rad/cm];|#Delta(dtbin/dr)| [tbin/cm]", m_connect_dphi_slope, m_connect_dtbin_slope).c_str(), 160, 0.0, std::max(0.08, 5.0 * m_connect_dphi_slope), 160, 0.0, std::max(20.0, 5.0 * m_connect_dtbin_slope));
0542   m_h_dphi_vs_dmphi = new TH2D("h_dphi_vs_dmphi", std::format("#phi position vs slope residual, cuts #Delta#phi<{:.4g} #Delta(d#phi/dr)<{:.4g};|#Delta#phi| [rad];|#Delta(d#phi/dr)| [rad/cm]", m_connect_dphi, m_connect_dphi_slope).c_str(), 160, 0.0, std::max(0.2, 5.0 * m_connect_dphi), 160, 0.0, std::max(0.08, 5.0 * m_connect_dphi_slope));
0543
0544   m_h_tbin_slope_vs_first_tbin = new TH2D("h_tbin_slope_vs_first_tbin", "tested connection dtbin/dr vs first timebin;first timebin;dtbin/dr [tbin/cm]", 200, 0.0, 600.0, 200, -20.0, 20.0);
0545   m_h_tbin_slope_vs_last_tbin = new TH2D("h_tbin_slope_vs_last_tbin", "tested connection dtbin/dr vs last timebin;last timebin;dtbin/dr [tbin/cm]", 200, 0.0, 600.0, 200, -20.0, 20.0);
0546   m_h_track_tbin_slope_vs_tbin_span_3modules = new TH2D("h_track_tbin_slope_vs_tbin_span_3modules", "3-module tracks dtbin/dr vs last-first timebin;last timebin - first timebin;dtbin/dr [tbin/cm]", 200, -600.0, 600.0, 200, -20.0, 20.0);
0547   m_h_track_tbin_slope_vs_first_tbin_3modules = new TH2D("h_track_tbin_slope_vs_first_tbin_3modules", "3-module tracks dtbin/dr vs first timebin;first timebin;dtbin/dr [tbin/cm]", 200, 0.0, 600.0, 200, -20.0, 20.0);
0548   m_h_track_tbin_slope_vs_last_tbin_3modules = new TH2D("h_track_tbin_slope_vs_last_tbin_3modules", "3-module tracks dtbin/dr vs last timebin;last timebin;dtbin/dr [tbin/cm]", 200, 0.0, 600.0, 200, -20.0, 20.0);
0549
0550   m_h_layer_gap = new TH1D("h_layer_gap", "accepted connection layer gap;b.first_layer - a.last_layer - 1;accepted connections", 16, -0.5, 15.5);
0551   m_h_nsegments = new TH1D("h_nsegments", "pieces per assembled track;nsegments;assembled tracks", 16, -0.5, 15.5);
0552   m_h_matched_sector_delta = new TH1D("h_matched_sector_delta", "accepted matched sector difference;wrapped #Delta sector;accepted connections", 25, -12.5, 12.5);
0553 }
0554
0555 void Tpc_AssembledTrackReco::write_debug_histograms()
0556 {
0557   if (!m_debugOutputFile)
0558   {
0559     return;
0560   }
0561   m_debugOutputFile->cd();
0562
0563   if (m_h_dphi)
0564   {
0565     m_h_dphi->Write();
0566   }
0567   if (m_h_dtbin)
0568   {
0569     m_h_dtbin->Write();
0570   }
0571   if (m_h_dmphi)
0572   {
0573     m_h_dmphi->Write();
0574   }
0575   if (m_h_dmtbin)
0576   {
0577     m_h_dmtbin->Write();
0578   }
0579   if (m_h_score)
0580   {
0581     m_h_score->Write();
0582   }
0583   if (m_h_dphi_vs_dtbin)
0584   {
0585     m_h_dphi_vs_dtbin->Write();
0586   }
0587   if (m_h_dmphi_vs_dmtbin)
0588   {
0589     m_h_dmphi_vs_dmtbin->Write();
0590   }
0591   if (m_h_dphi_vs_dmphi)
0592   {
0593     m_h_dphi_vs_dmphi->Write();
0594   }
0595   if (m_h_tbin_slope_vs_first_tbin)
0596   {
0597     m_h_tbin_slope_vs_first_tbin->Write();
0598   }
0599   if (m_h_tbin_slope_vs_last_tbin)
0600   {
0601     m_h_tbin_slope_vs_last_tbin->Write();
0602   }
0603   if (m_h_track_tbin_slope_vs_tbin_span_3modules)
0604   {
0605     m_h_track_tbin_slope_vs_tbin_span_3modules->Write();
0606   }
0607   if (m_h_track_tbin_slope_vs_first_tbin_3modules)
0608   {
0609     m_h_track_tbin_slope_vs_first_tbin_3modules->Write();
0610   }
0611   if (m_h_track_tbin_slope_vs_last_tbin_3modules)
0612   {
0613     m_h_track_tbin_slope_vs_last_tbin_3modules->Write();
0614   }
0615   if (m_h_layer_gap)
0616   {
0617     m_h_layer_gap->Write();
0618   }
0619   if (m_h_nsegments)
0620   {
0621     m_h_nsegments->Write();
0622   }
0623   if (m_h_matched_sector_delta)
0624   {
0625     m_h_matched_sector_delta->Write();
0626   }
0627 }
0628
0629 void Tpc_AssembledTrackReco::reset_tree_vars()
0630 {
0631   m_tree_event = m_event;
0632   m_tree_track_id.clear();
0633   m_tree_side.clear();
0634   m_tree_nsegments.clear();
0635   m_tree_nblobs.clear();
0636   m_tree_nrawhits.clear();
0637   m_tree_first_layer.clear();
0638   m_tree_last_layer.clear();
0639   m_tree_first_sector.clear();
0640   m_tree_last_sector.clear();
0641   m_tree_first_region.clear();
0642   m_tree_last_region.clear();
0643   m_tree_source_assembled_track_id.clear();
0644   m_tree_source_inmodule_track_id.clear();
0645   m_tree_source_region.clear();
0646   m_tree_source_sector.clear();
0647   m_tree_source_side.clear();
0648   m_tree_hit_assembled_track_id.clear();
0649   m_tree_hit_hitsetkey.clear();
0650   m_tree_hit_hitkey.clear();
0651 }
0652
0653 bool Tpc_AssembledTrackReco::make_piece(unsigned int source_index, Piece& p) const
0654 {
0655   const Tpc_ModuleTrack* trk = m_tpcModuleTrackContainer->get_track(source_index);
0656   if (!trk || !trk->isValid())
0657   {
0658     return false;
0659   }
0660   if (trk->get_last_layer() < trk->get_first_layer())
0661   {
0662     return false;
0663   }
0664   if (!m_idealPadMap || !m_idealPadMap->is_loaded() || !m_hits)
0665   {
0666     return false;
0667   }
0668
0669   p.source_index = source_index;
0670   p.source_track_id = trk->get_track_id();
0671   p.event = trk->get_event();
0672   p.region = trk->get_region();
0673   p.sector = trk->get_sector();
0674   p.side = trk->get_side();
0675   p.first_layer = trk->get_first_layer();
0676   p.last_layer = trk->get_last_layer();
0677   p.nblobs = trk->get_nblobs();
0678   p.nrawhits = trk->get_nrawhits();
0679
0680   p.radius_values.clear();
0681   p.phi_values.clear();
0682   p.tbin_values.clear();
0683   p.weights.clear();
0684   p.hitsetkeys.clear();
0685   p.hitkeys.clear();
0686
0687   double maxadc = 0.0;
0688   for (unsigned int ih = 0; ih < trk->size_hit_indices(); ++ih)
0689   {
0690     const Tpc_ModuleTrack::HitIndex hi = trk->get_hit_index(ih);
0691     TrkrHitSet* hitset = m_hits->findHitSet(hi.first);
0692     if (!hitset)
0693     {
0694       continue;
0695     }
0696     TrkrHit* hit = hitset->getHit(hi.second);
0697     if (!hit)
0698     {
0699       continue;
0700     }
0701     const double adc = static_cast<double>(hit->getAdc());
0702     maxadc = std::max(adc, maxadc);
0703   }
0704
0705   for (unsigned int ih = 0; ih < trk->size_hit_indices(); ++ih)
0706   {
0707     const Tpc_ModuleTrack::HitIndex hi = trk->get_hit_index(ih);
0708     TrkrHitSet* hitset = m_hits->findHitSet(hi.first);
0709     if (!hitset)
0710     {
0711       continue;
0712     }
0713     TrkrHit* hit = hitset->getHit(hi.second);
0714     if (!hit)
0715     {
0716       continue;
0717     }
0718
0719     const unsigned int layer = TrkrDefs::getLayer(hi.first);
0720     const unsigned int pad = TpcDefs::getPad(hi.second);
0721     const unsigned int tbin = TpcDefs::getTBin(hi.second);
0722
0723     const double radius = m_idealPadMap->get_radius(layer);
0724     const double phi = wrap_to_pi(m_idealPadMap->get_phi(static_cast<unsigned int>(p.side), layer, pad));
0725     if (!std::isfinite(radius) || !std::isfinite(phi))
0726     {
0727       continue;
0728     }
0729
0730     p.radius_values.push_back(radius);
0731     p.phi_values.push_back(phi);
0732     p.tbin_values.push_back(static_cast<double>(tbin));
0733     p.weights.push_back(Tpc_FittingTools::adcWeight(static_cast<double>(hit->getAdc()), maxadc, 0.5, 0.15));
0734     p.hitsetkeys.push_back(hi.first);
0735     p.hitkeys.push_back(hi.second);
0736   }
0737
0738   if (p.radius_values.size() < 2)
0739   {
0740     return false;
0741   }
0742
0743   std::vector<unsigned int> order;
0744   order.reserve(p.radius_values.size());
0745   for (unsigned int i = 0; i < p.radius_values.size(); ++i)
0746   {
0747     order.push_back(i);
0748   }
0749   std::sort(order.begin(), order.end(), RadiusSort(&p.radius_values));
0750
0751   std::vector<double> r_sorted;
0752   std::vector<double> phi_sorted;
0753   std::vector<double> tbin_sorted;
0754   std::vector<double> w_sorted;
0755   std::vector<TrkrDefs::hitsetkey> hsk_sorted;
0756   std::vector<TrkrDefs::hitkey> hk_sorted;
0757   r_sorted.reserve(order.size());
0758   phi_sorted.reserve(order.size());
0759   tbin_sorted.reserve(order.size());
0760   w_sorted.reserve(order.size());
0761   hsk_sorted.reserve(order.size());
0762   hk_sorted.reserve(order.size());
0763
0764   for (unsigned int i : order)
0765   {
0766     double phi = p.phi_values[i];
0767     if (!phi_sorted.empty())
0768     {
0769       phi = unwrap_phi_to_reference(phi, phi_sorted.back());
0770     }
0771     r_sorted.push_back(p.radius_values[i]);
0772     phi_sorted.push_back(phi);
0773     tbin_sorted.push_back(p.tbin_values[i]);
0774     w_sorted.push_back(p.weights[i]);
0775     hsk_sorted.push_back(p.hitsetkeys[i]);
0776     hk_sorted.push_back(p.hitkeys[i]);
0777   }
0778
0779   p.radius_values.swap(r_sorted);
0780   p.phi_values.swap(phi_sorted);
0781   p.tbin_values.swap(tbin_sorted);
0782   p.weights.swap(w_sorted);
0783   p.hitsetkeys.swap(hsk_sorted);
0784   p.hitkeys.swap(hk_sorted);
0785
0786   double chi2_phi = 0.0;
0787   double chi2_tbin = 0.0;
0788   int ndof_phi = 0;
0789   int ndof_tbin = 0;
0790   return fit_points(p.radius_values, p.phi_values, p.tbin_values, p.weights, m_useSagittaPhiFit,
0791                     p.phi_slope, p.phi_intercept, p.phi_S, p.phi_x0, p.phi_invR,
0792                     p.phi_theta, p.phi_bline, p.phi_sagitta_ok,
0793                     p.tbin_slope, p.tbin_intercept, chi2_phi, chi2_tbin, ndof_phi, ndof_tbin);
0794 }
0795
0796 double Tpc_AssembledTrackReco::predict_phi(const Piece& p, double radius) const
0797 {
0798   if (m_useSagittaPhiFit && p.phi_sagitta_ok)
0799   {
0800     return predict_sagitta_phi(radius, p.phi_S, p.phi_x0, p.phi_invR, p.phi_theta, p.phi_bline);
0801   }
0802   return p.phi_slope * radius + p.phi_intercept;
0803 }
0804
0805 double Tpc_AssembledTrackReco::predict_phi(const Candidate& c, double radius) const
0806 {
0807   if (m_useSagittaPhiFit && c.phi_sagitta_ok)
0808   {
0809     return predict_sagitta_phi(radius, c.phi_S, c.phi_x0, c.phi_invR, c.phi_theta, c.phi_bline);
0810   }
0811   return c.phi_slope * radius + c.phi_intercept;
0812 }
0813
0814 double Tpc_AssembledTrackReco::predict_phi_slope(const Piece& p, double radius) const
0815 {
0816   if (!(m_useSagittaPhiFit && p.phi_sagitta_ok))
0817   {
0818     return p.phi_slope;
0819   }
0820   const double eps = 1.0e-3;
0821   return (predict_phi(p, radius + eps) - predict_phi(p, radius - eps)) / (2.0 * eps);
0822 }
0823
0824 double Tpc_AssembledTrackReco::predict_phi_slope(const Candidate& c, double radius) const
0825 {
0826   if (!(m_useSagittaPhiFit && c.phi_sagitta_ok))
0827   {
0828     return c.phi_slope;
0829   }
0830   const double eps = 1.0e-3;
0831   return (predict_phi(c, radius + eps) - predict_phi(c, radius - eps)) / (2.0 * eps);
0832 }
0833
0834 bool Tpc_AssembledTrackReco::refit_candidate(const std::vector<Piece>& pieces, const std::vector<unsigned int>& piece_indices, Candidate& c) const
0835 {
0836   if (piece_indices.empty())
0837   {
0838     return false;
0839   }
0840
0841   std::vector<double> radius;
0842   std::vector<double> phi;
0843   std::vector<double> tbin;
0844   std::vector<double> weight;
0845   c = Candidate();
0846   c.piece_indices = piece_indices;
0847
0848   for (unsigned int ii = 0; ii < piece_indices.size(); ++ii)
0849   {
0850     const Piece& p = pieces[piece_indices[ii]];
0851
0852     for (unsigned int ih = 0; ih < p.radius_values.size(); ++ih)
0853     {
0854       double phiv = p.phi_values[ih];
0855       if (!phi.empty())
0856       {
0857         phiv = unwrap_phi_to_reference(phiv, phi.back());
0858       }
0859       radius.push_back(p.radius_values[ih]);
0860       phi.push_back(phiv);
0861       tbin.push_back(p.tbin_values[ih]);
0862       weight.push_back(p.weights[ih]);
0863     }
0864
0865     if (ii == 0)
0866     {
0867       c.event = p.event;
0868       c.side = p.side;
0869       c.first_layer = p.first_layer;
0870       c.last_layer = p.last_layer;
0871       c.first_sector = p.sector;
0872       c.last_sector = p.sector;
0873       c.first_region = p.region;
0874       c.last_region = p.region;
0875     }
0876     else
0877     {
0878       if (p.first_layer < c.first_layer)
0879       {
0880         c.first_layer = p.first_layer;
0881         c.first_sector = p.sector;
0882         c.first_region = p.region;
0883       }
0884       if (p.last_layer > c.last_layer)
0885       {
0886         c.last_layer = p.last_layer;
0887         c.last_sector = p.sector;
0888         c.last_region = p.region;
0889       }
0890     }
0891
0892     c.nblobs += p.nblobs;
0893     c.nrawhits += p.nrawhits;
0894     for (unsigned int ih = 0; ih < p.hitsetkeys.size(); ++ih)
0895     {
0896       c.hitsetkeys.push_back(p.hitsetkeys[ih]);
0897       c.hitkeys.push_back(p.hitkeys[ih]);
0898     }
0899   }
0900
0901   c.nsegments = static_cast<unsigned int>(piece_indices.size());
0902   return fit_points(radius, phi, tbin, weight, m_useSagittaPhiFit,
0903                     c.phi_slope, c.phi_intercept, c.phi_S, c.phi_x0, c.phi_invR,
0904                     c.phi_theta, c.phi_bline, c.phi_sagitta_ok,
0905                     c.tbin_slope_r, c.tbin_intercept_r,
0906                     c.chi2_phi, c.chi2_tbin, c.ndof_phi, c.ndof_tbin);
0907 }
0908
0909 bool Tpc_AssembledTrackReco::candidates_can_connect(const Candidate& a, const Piece& b, double& score, double& b_phi_intercept_shifted) const
0910 {
0911   score = std::numeric_limits<double>::max();
0912   b_phi_intercept_shifted = b.phi_intercept;
0913
0914   if (a.side != b.side)
0915   {
0916     return false;
0917   }
0918   if (a.last_layer >= b.first_layer)
0919   {
0920     return false;
0921   }
0922
0923   const unsigned int gap = b.first_layer - a.last_layer - 1;
0924   if (gap > m_connectMaxLayerGap)
0925   {
0926     return false;
0927   }
0928
0929   const double ra = m_idealPadMap->get_radius(a.last_layer);
0930   const double rb = m_idealPadMap->get_radius(b.first_layer);
0931   if (!std::isfinite(ra) || !std::isfinite(rb) || ra <= 0.0 || rb <= 0.0)
0932   {
0933     return false;
0934   }
0935
0936   const double rmatch = 0.5 * (ra + rb);
0937   const double phi_a = predict_phi(a, rmatch);
0938   const double phi_b_raw = predict_phi(b, rmatch);
0939   const double phi_b = unwrap_phi_to_reference(phi_b_raw, phi_a);  // NOLINT(readability-suspicious-call-argument)
0940   b_phi_intercept_shifted = b.phi_intercept + (phi_b - phi_b_raw);
0941
0942   const double tbin_a = a.tbin_slope_r * rmatch + a.tbin_intercept_r;
0943   const double tbin_b = b.tbin_slope * rmatch + b.tbin_intercept;
0944   const double dphi = std::fabs(phi_a - phi_b);
0945   const double dtbin = std::fabs(tbin_a - tbin_b);
0946   const double dmphi = std::fabs(predict_phi_slope(a, rmatch) - predict_phi_slope(b, rmatch));
0947   const double dmtbin = std::fabs(a.tbin_slope_r - b.tbin_slope);
0948
0949   {
0950     std::lock_guard<std::mutex> lock(m_debugMutex);
0951     if (m_h_dphi)
0952     {
0953       m_h_dphi->Fill(dphi);
0954     }
0955     if (m_h_dtbin)
0956     {
0957       m_h_dtbin->Fill(dtbin);
0958     }
0959     if (m_h_dmphi)
0960     {
0961       m_h_dmphi->Fill(dmphi);
0962     }
0963     if (m_h_dmtbin)
0964     {
0965       m_h_dmtbin->Fill(dmtbin);
0966     }
0967     if (m_h_dphi_vs_dtbin)
0968     {
0969       m_h_dphi_vs_dtbin->Fill(dphi, dtbin);
0970     }
0971     if (m_h_dmphi_vs_dmtbin)
0972     {
0973       m_h_dmphi_vs_dmtbin->Fill(dmphi, dmtbin);
0974     }
0975     if (m_h_dphi_vs_dmphi)
0976     {
0977       m_h_dphi_vs_dmphi->Fill(dphi, dmphi);
0978     }
0979     if (m_h_tbin_slope_vs_last_tbin)
0980     {
0981       m_h_tbin_slope_vs_last_tbin->Fill(a.tbin_slope_r * ra + a.tbin_intercept_r, a.tbin_slope_r);
0982     }
0983     if (m_h_tbin_slope_vs_first_tbin && !b.tbin_values.empty())
0984     {
0985       m_h_tbin_slope_vs_first_tbin->Fill(b.tbin_values.front(), b.tbin_slope);
0986     }
0987   }
0988
0989   if (dphi > m_connect_dphi)
0990   {
0991     return false;
0992   }
0993   if (dtbin > m_connect_dtbin)
0994   {
0995     return false;
0996   }
0997   if (dmphi > m_connect_dphi_slope)
0998   {
0999     return false;
1000   }
1001   if (dmtbin > m_connect_dtbin_slope)
1002   {
1003     return false;
1004   }
1005
1006   constexpr double w_phi = 1.0;    // φ position — highest weight
1007   constexpr double w_mphi = 1.0;   // φ slope
1008   constexpr double w_tbin = 1.0;   // tbin position — softer penalty
1009   constexpr double w_mtbin = 2.0;  // tbin slope
1010
1011   score = w_phi * (dphi / m_connect_dphi) * (dphi / m_connect_dphi) + w_tbin * (dtbin / m_connect_dtbin) * (dtbin / m_connect_dtbin) + w_mphi * (dmphi / m_connect_dphi_slope) * (dmphi / m_connect_dphi_slope) + w_mtbin * (dmtbin / m_connect_dtbin_slope) * (dmtbin / m_connect_dtbin_slope) + 0.05 * static_cast<double>(gap);
1012
1013   return true;
1014 }
1015
1016 void Tpc_AssembledTrackReco::connect_sector_pieces(const std::vector<Piece>& pieces, int side, unsigned int sector, std::vector<Candidate>& output) const
1017 {
1018   std::vector<unsigned int> order;
1019   for (unsigned int i = 0; i < pieces.size(); ++i)
1020   {
1021     if (pieces[i].side == side && pieces[i].sector == sector)
1022     {
1023       order.push_back(i);
1024     }
1025   }
1026   if (order.empty())
1027   {
1028     return;
1029   }
1030
1031   std::sort(order.begin(), order.end(), PieceStartSort(&pieces));
1032   std::vector<int> used(pieces.size(), 0);
1033
1034   for (unsigned int io = 0; io < order.size(); ++io)
1035   {
1036     const unsigned int iseed = order[io];
1037     if (used[iseed])
1038     {
1039       continue;
1040     }
1041
1042     std::vector<unsigned int> current_indices;
1043     current_indices.push_back(iseed);
1044     used[iseed] = 1;
1045
1046     Candidate current;
1047     if (!refit_candidate(pieces, current_indices, current))
1048     {
1049       continue;
1050     }
1051
1052     bool merged_any = true;
1053     while (merged_any)
1054     {
1055       merged_any = false;
1056       int best_j = -1;
1057       double best_score = std::numeric_limits<double>::max();
1058
1059       for (unsigned int j : order)
1060       {
1061         if (used[j])
1062         {
1063           continue;
1064         }
1065
1066         double score = 0.0;
1067         double shifted_intercept = 0.0;
1068         if (!candidates_can_connect(current, pieces[j], score, shifted_intercept))
1069         {
1070           continue;
1071         }
1072
1073         if (score < best_score)
1074         {
1075           best_score = score;
1076           best_j = static_cast<int>(j);
1077         }
1078       }
1079
1080       if (best_j >= 0)
1081       {
1082         std::vector<unsigned int> trial_indices = current_indices;
1083         trial_indices.push_back(static_cast<unsigned int>(best_j));
1084
1085         Candidate refit;
1086         if (refit_candidate(pieces, trial_indices, refit))
1087         {
1088           const Piece& accepted_piece = pieces[static_cast<unsigned int>(best_j)];
1089           const unsigned int accepted_gap = accepted_piece.first_layer - current.last_layer - 1;
1090           {
1091             std::lock_guard<std::mutex> lock(m_debugMutex);
1092             if (m_h_score)
1093             {
1094               m_h_score->Fill(best_score);
1095             }
1096             if (m_h_layer_gap)
1097             {
1098               m_h_layer_gap->Fill(static_cast<double>(accepted_gap));
1099             }
1100             if (m_h_matched_sector_delta)
1101             {
1102               m_h_matched_sector_delta->Fill(static_cast<double>(wrapped_sector_delta(current.last_sector, accepted_piece.sector)));
1103             }
1104           }
1105
1106           current = refit;
1107           current_indices.swap(trial_indices);
1108           used[best_j] = 1;
1109           merged_any = true;
1110         }
1111       }
1112     }
1113
1114     output.push_back(current);
1115   }
1116 }
1117
1118 bool Tpc_AssembledTrackReco::candidates_can_connect(const Candidate& a, const Candidate& b, double& score) const
1119 {
1120   score = std::numeric_limits<double>::max();
1121
1122   if (a.side != b.side)
1123   {
1124     return false;
1125   }
1126   if (a.last_layer >= b.first_layer)
1127   {
1128     return false;
1129   }
1130   if (a.last_sector == b.first_sector)
1131   {
1132     return false;
1133   }
1134
1135   const unsigned int gap = b.first_layer - a.last_layer - 1;
1136   if (gap > m_connectMaxLayerGap)
1137   {
1138     return false;
1139   }
1140
1141   const double ra = m_idealPadMap->get_radius(a.last_layer);
1142   const double rb = m_idealPadMap->get_radius(b.first_layer);
1143   if (!std::isfinite(ra) || !std::isfinite(rb) || ra <= 0.0 || rb <= 0.0)
1144   {
1145     return false;
1146   }
1147
1148   const double rmatch = 0.5 * (ra + rb);
1149   const double phi_a = predict_phi(a, rmatch);
1150   const double phi_b = unwrap_phi_to_reference(predict_phi(b, rmatch), phi_a);
1151   const double tbin_a = a.tbin_slope_r * rmatch + a.tbin_intercept_r;
1152   const double tbin_b = b.tbin_slope_r * rmatch + b.tbin_intercept_r;
1153   const double dphi = std::fabs(phi_a - phi_b);
1154   const double dtbin = std::fabs(tbin_a - tbin_b);
1155   const double dmphi = std::fabs(predict_phi_slope(a, rmatch) - predict_phi_slope(b, rmatch));
1156   const double dmtbin = std::fabs(a.tbin_slope_r - b.tbin_slope_r);
1157
1158   {
1159     std::lock_guard<std::mutex> lock(m_debugMutex);
1160     if (m_h_dphi)
1161     {
1162       m_h_dphi->Fill(dphi);
1163     }
1164     if (m_h_dtbin)
1165     {
1166       m_h_dtbin->Fill(dtbin);
1167     }
1168     if (m_h_dmphi)
1169     {
1170       m_h_dmphi->Fill(dmphi);
1171     }
1172     if (m_h_dmtbin)
1173     {
1174       m_h_dmtbin->Fill(dmtbin);
1175     }
1176     if (m_h_dphi_vs_dtbin)
1177     {
1178       m_h_dphi_vs_dtbin->Fill(dphi, dtbin);
1179     }
1180     if (m_h_dmphi_vs_dmtbin)
1181     {
1182       m_h_dmphi_vs_dmtbin->Fill(dmphi, dmtbin);
1183     }
1184     if (m_h_dphi_vs_dmphi)
1185     {
1186       m_h_dphi_vs_dmphi->Fill(dphi, dmphi);
1187     }
1188     if (m_h_tbin_slope_vs_last_tbin)
1189     {
1190       m_h_tbin_slope_vs_last_tbin->Fill(a.tbin_slope_r * ra + a.tbin_intercept_r, a.tbin_slope_r);
1191     }
1192     if (m_h_tbin_slope_vs_first_tbin)
1193     {
1194       m_h_tbin_slope_vs_first_tbin->Fill(b.tbin_slope_r * rb + b.tbin_intercept_r, b.tbin_slope_r);
1195     }
1196   }
1197
1198   if (dphi > m_connect_dphi)
1199   {
1200     return false;
1201   }
1202   if (dtbin > m_connect_dtbin)
1203   {
1204     return false;
1205   }
1206   if (dmphi > m_connect_dphi_slope)
1207   {
1208     return false;
1209   }
1210   if (dmtbin > m_connect_dtbin_slope)
1211   {
1212     return false;
1213   }
1214
1215   constexpr double w_phi = 1.0;
1216   constexpr double w_mphi = 1.0;
1217   constexpr double w_tbin = 1.0;
1218   constexpr double w_mtbin = 2.0;
1219
1220   score = w_phi * (dphi / m_connect_dphi) * (dphi / m_connect_dphi) + w_tbin * (dtbin / m_connect_dtbin) * (dtbin / m_connect_dtbin) + w_mphi * (dmphi / m_connect_dphi_slope) * (dmphi / m_connect_dphi_slope) + w_mtbin * (dmtbin / m_connect_dtbin_slope) * (dmtbin / m_connect_dtbin_slope) + 0.05 * static_cast<double>(gap);
1221
1222   return true;
1223 }
1224
1225 void Tpc_AssembledTrackReco::connect_side_candidates(const std::vector<Piece>& pieces,
1226                                                      const std::vector<Candidate>& seeds,
1227                                                      int side,
1228                                                      std::vector<Candidate>& output) const
1229 {
1230   std::vector<unsigned int> order;
1231   for (unsigned int i = 0; i < seeds.size(); ++i)
1232   {
1233     if (seeds[i].side == side)
1234     {
1235       order.push_back(i);
1236     }
1237   }
1238   if (order.empty())
1239   {
1240     return;
1241   }
1242
1243   std::sort(order.begin(), order.end(), CandidateStartSort(&seeds));
1244   std::vector<int> used(seeds.size(), 0);
1245
1246   for (unsigned int io = 0; io < order.size(); ++io)
1247   {
1248     const unsigned int iseed = order[io];
1249     if (used[iseed])
1250     {
1251       continue;
1252     }
1253
1254     std::vector<unsigned int> current_indices = seeds[iseed].piece_indices;
1255     used[iseed] = 1;
1256
1257     Candidate current;
1258     if (!refit_candidate(pieces, current_indices, current))
1259     {
1260       continue;
1261     }
1262
1263     bool merged_any = true;
1264     while (merged_any)
1265     {
1266       merged_any = false;
1267       int best_j = -1;
1268       double best_score = std::numeric_limits<double>::max();
1269
1270       for (unsigned int j : order)
1271       {
1272         if (used[j])
1273         {
1274           continue;
1275         }
1276
1277         double score = 0.0;
1278         if (!candidates_can_connect(current, seeds[j], score))
1279         {
1280           continue;
1281         }
1282
1283         if (score < best_score)
1284         {
1285           best_score = score;
1286           best_j = static_cast<int>(j);
1287         }
1288       }
1289
1290       if (best_j >= 0)
1291       {
1292         std::vector<unsigned int> trial_indices = current_indices;
1293         const Candidate& accepted_seed = seeds[static_cast<unsigned int>(best_j)];
1294         trial_indices.insert(trial_indices.end(), accepted_seed.piece_indices.begin(), accepted_seed.piece_indices.end());
1295
1296         Candidate refit;
1297         if (refit_candidate(pieces, trial_indices, refit))
1298         {
1299           const unsigned int accepted_gap = accepted_seed.first_layer - current.last_layer - 1;
1300           {
1301             std::lock_guard<std::mutex> lock(m_debugMutex);
1302             if (m_h_score)
1303             {
1304               m_h_score->Fill(best_score);
1305             }
1306             if (m_h_layer_gap)
1307             {
1308               m_h_layer_gap->Fill(static_cast<double>(accepted_gap));
1309             }
1310             if (m_h_matched_sector_delta)
1311             {
1312               m_h_matched_sector_delta->Fill(static_cast<double>(wrapped_sector_delta(current.last_sector, accepted_seed.first_sector)));
1313             }
1314           }
1315
1316           current = refit;
1317           current_indices.swap(trial_indices);
1318           used[best_j] = 1;
1319           merged_any = true;
1320         }
1321       }
1322     }
1323
1324     output.push_back(current);
1325   }
1326 }
1327
1328 Tpc_AssembledTrackReco::SeedParameters
1329 Tpc_AssembledTrackReco::make_seed_parameters(const Candidate& c) const
1330 {
1331   SeedParameters seed;
1332   if (!m_idealPadMap)
1333   {
1334     return seed;
1335   }
1336
1337   const double r_first = m_idealPadMap->get_radius(c.first_layer);
1338   const double r_last = m_idealPadMap->get_radius(c.last_layer);
1339   if (!std::isfinite(r_first) || !std::isfinite(r_last))
1340   {
1341     return seed;
1342   }
1343
1344   const double radius = 0.5 * (r_first + r_last);
1345   const double phi = wrap_to_pi(predict_phi(c, radius));
1346   const double dphi_dr = predict_phi_slope(c, radius);
1347   const double dtbin_dr = c.tbin_slope_r;
1348   const double tbin = c.tbin_slope_r * radius + c.tbin_intercept_r;
1349   if (!std::isfinite(radius) || !std::isfinite(phi) ||
1350       !std::isfinite(dphi_dr) || !std::isfinite(dtbin_dr) || !std::isfinite(tbin))
1351   {
1352     return seed;
1353   }
1354
1355   seed.x = radius * std::cos(phi);
1356   seed.y = radius * std::sin(phi);
1357   seed.z = tbin;
1358
1359   const double dx_dr = std::cos(phi) - radius * std::sin(phi) * dphi_dr;
1360   const double dy_dr = std::sin(phi) + radius * std::cos(phi) * dphi_dr;
1361   const double dz_dr = dtbin_dr;
1362   const double norm = std::sqrt(dx_dr * dx_dr + dy_dr * dy_dr + dz_dr * dz_dr);
1363   if (!std::isfinite(norm) || norm <= 0.0)
1364   {
1365     return seed;
1366   }
1367
1368   seed.px = dx_dr / norm;
1369   seed.py = dy_dr / norm;
1370   seed.pz = dz_dr / norm;
1371
1372   for (auto& i : seed.cov)
1373   {
1374     for (double& j : i)
1375     {
1376       j = 0.0;
1377     }
1378   }
1379   seed.cov[0][0] = m_seedSigmaX * m_seedSigmaX;
1380   seed.cov[1][1] = m_seedSigmaY * m_seedSigmaY;
1381   seed.cov[2][2] = m_seedSigmaZ * m_seedSigmaZ;
1382   seed.cov[3][3] = m_seedSigmaPx * m_seedSigmaPx;
1383   seed.cov[4][4] = m_seedSigmaPy * m_seedSigmaPy;
1384   seed.cov[5][5] = m_seedSigmaPz * m_seedSigmaPz;
1385   seed.ok = std::isfinite(seed.x) && std::isfinite(seed.y) && std::isfinite(seed.z) &&
1386             std::isfinite(seed.px) && std::isfinite(seed.py) && std::isfinite(seed.pz);
1387   return seed;
1388 }
1389
1390 int Tpc_AssembledTrackReco::process_event(PHCompositeNode* /*unused*/)
1391 {
1392   reset_tree_vars();
1393   if (m_assembledTrackContainer)
1394   {
1395     m_assembledTrackContainer->Reset();
1396   }
1397
1398   std::vector<Piece> pieces;
1399   if (m_tpcModuleTrackContainer)
1400   {
1401     const unsigned int n = m_tpcModuleTrackContainer->size();
1402     pieces.reserve(n);
1403     for (unsigned int i = 0; i < n; ++i)
1404     {
1405       Piece p;
1406       if (make_piece(i, p))
1407       {
1408         pieces.push_back(p);
1409       }
1410     }
1411   }
1412
1413   std::vector<std::vector<Candidate> > sector_outputs(24);
1414   std::vector<std::thread> workers;
1415   workers.reserve(24);
1416   for (int side = 0; side < 2; ++side)
1417   {
1418     for (unsigned int sector = 0; sector < 12; ++sector)
1419     {
1420       const unsigned int index = static_cast<unsigned int>(side) * 12 + sector;
1421       workers.emplace_back(&Tpc_AssembledTrackReco::connect_sector_pieces, this, std::cref(pieces), side, sector, std::ref(sector_outputs[index]));
1422     }
1423   }
1424   for (std::thread& worker : workers)
1425   {
1426     worker.join();
1427   }
1428
1429   std::vector<Candidate> sector_tracks;
1430   for (auto& sector_output : sector_outputs)
1431   {
1432     sector_tracks.insert(sector_tracks.end(), sector_output.begin(), sector_output.end());
1433   }
1434
1435   std::vector<Candidate> assembled_tracks;
1436   connect_side_candidates(pieces, sector_tracks, 0, assembled_tracks);
1437   connect_side_candidates(pieces, sector_tracks, 1, assembled_tracks);
1438
1439   for (unsigned int it = 0; it < assembled_tracks.size(); ++it)
1440   {
1441     const Candidate& c = assembled_tracks[it];
1442     const unsigned int assembled_id = m_assembledTrackContainer ? m_assembledTrackContainer->size() : it;
1443     if (m_h_nsegments)
1444     {
1445       m_h_nsegments->Fill(static_cast<double>(c.nsegments));
1446     }
1447     if (c.nsegments == 3 && m_idealPadMap)
1448     {
1449       const double rfirst = m_idealPadMap->get_radius(c.first_layer);
1450       const double rlast = m_idealPadMap->get_radius(c.last_layer);
1451       if (std::isfinite(rfirst) && std::isfinite(rlast))
1452       {
1453         const double first_tbin = c.tbin_slope_r * rfirst + c.tbin_intercept_r;
1454         const double last_tbin = c.tbin_slope_r * rlast + c.tbin_intercept_r;
1455         const double dtbin_track = last_tbin - first_tbin;
1456         if (m_h_track_tbin_slope_vs_tbin_span_3modules)
1457         {
1458           m_h_track_tbin_slope_vs_tbin_span_3modules->Fill(dtbin_track, c.tbin_slope_r);
1459         }
1460         if (m_h_track_tbin_slope_vs_first_tbin_3modules)
1461         {
1462           m_h_track_tbin_slope_vs_first_tbin_3modules->Fill(first_tbin, c.tbin_slope_r);
1463         }
1464         if (m_h_track_tbin_slope_vs_last_tbin_3modules)
1465         {
1466           m_h_track_tbin_slope_vs_last_tbin_3modules->Fill(last_tbin, c.tbin_slope_r);
1467         }
1468       }
1469     }
1470
1471     Tpc_AssembledTrackv1* out = new Tpc_AssembledTrackv1();
1472     out->set_event(static_cast<unsigned int>(m_event));
1473     out->set_track_id(assembled_id);
1474     out->set_side(c.side);
1475     out->set_nsegments(c.nsegments);
1476     out->set_nblobs(c.nblobs);
1477     out->set_nrawhits(c.nrawhits);
1478     out->set_first_layer(c.first_layer);
1479     out->set_last_layer(c.last_layer);
1480     out->set_first_sector(c.first_sector);
1481     out->set_last_sector(c.last_sector);
1482     out->set_first_region(c.first_region);
1483     out->set_last_region(c.last_region);
1484
1485     const SeedParameters seed = make_seed_parameters(c);
1486     if (seed.ok)
1487     {
1488       out->set_seed_valid(1);
1489       out->set_seed_x(seed.x);
1490       out->set_seed_y(seed.y);
1491       out->set_seed_z(seed.z);
1492       out->set_seed_px(seed.px);
1493       out->set_seed_py(seed.py);
1494       out->set_seed_pz(seed.pz);
1495       for (unsigned int iseed = 0; iseed < 6; ++iseed)
1496       {
1497         for (unsigned int jseed = 0; jseed < 6; ++jseed)
1498         {
1499           out->set_seed_cov(iseed, jseed, seed.cov[iseed][jseed]);
1500         }
1501       }
1502     }
1503
1504     m_tree_track_id.push_back(assembled_id);
1505     m_tree_side.push_back(c.side);
1506     m_tree_nsegments.push_back(c.nsegments);
1507     m_tree_nblobs.push_back(c.nblobs);
1508     m_tree_nrawhits.push_back(c.nrawhits);
1509     m_tree_first_layer.push_back(c.first_layer);
1510     m_tree_last_layer.push_back(c.last_layer);
1511     m_tree_first_sector.push_back(c.first_sector);
1512     m_tree_last_sector.push_back(c.last_sector);
1513     m_tree_first_region.push_back(c.first_region);
1514     m_tree_last_region.push_back(c.last_region);
1515
1516     for (unsigned int piece_indice : c.piece_indices)
1517     {
1518       const Piece& p = pieces[piece_indice];
1519       out->add_source_track(p.source_track_id, p.region, p.sector);
1520       m_tree_source_assembled_track_id.push_back(assembled_id);
1521       m_tree_source_inmodule_track_id.push_back(p.source_track_id);
1522       m_tree_source_region.push_back(p.region);
1523       m_tree_source_sector.push_back(p.sector);
1524       m_tree_source_side.push_back(p.side);
1525     }
1526
1527     for (unsigned int ih = 0; ih < c.hitsetkeys.size(); ++ih)
1528     {
1529       out->add_hit_index(c.hitsetkeys[ih], c.hitkeys[ih]);
1530       m_tree_hit_assembled_track_id.push_back(assembled_id);
1531       m_tree_hit_hitsetkey.push_back(static_cast<unsigned long long>(c.hitsetkeys[ih]));
1532       m_tree_hit_hitkey.push_back(static_cast<unsigned long long>(c.hitkeys[ih]));
1533     }
1534
1535     if (m_assembledTrackContainer)
1536     {
1537       m_assembledTrackContainer->add_track(out);
1538     }
1539     else
1540     {
1541       delete out;
1542     }
1543   }
1544
1545   if (m_tree)
1546   {
1547     m_tree->Fill();
1548   }
1549
1550   if (Verbosity() > 0)
1551   {
1552     std::cout << Name() << "::process_event - event " << m_event
1553               << " input pieces=" << pieces.size()
1554               << " assembled_tracks=" << assembled_tracks.size() << std::endl;
1555   }
1556
1557   ++m_event;
1558   return Fun4AllReturnCodes::EVENT_OK;
1559 }