File indexing completed on 2026-08-31 08:21:25
0001 #include "TpcCrossingFinder.h"
0002
0003 #include "IdealPadMap.h"
0004 #include "Tpc_AssembledTrack.h"
0005 #include "Tpc_AssembledTrackContainer.h"
0006 #include "TpcCrossingDecisionContainerv1.h"
0007 #include "TpcCrossingDecisionv1.h"
0008 #include "Tpc_FittingTools.h"
0009
0010 #include <cdbobjects/CDBTTree.h>
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/PHIODataNode.h>
0015 #include <phool/PHNodeIterator.h>
0016 #include <phool/PHObject.h>
0017 #include <phool/getClass.h>
0018
0019 #include <ffamodules/CDBInterface.h>
0020
0021 #include <globalvertex/SvtxVertex.h>
0022 #include <globalvertex/SvtxVertexMap.h>
0023 #include <trackbase/InttDefs.h>
0024 #include <trackbase/TpcDefs.h>
0025 #include <trackbase/TrkrClusterContainer.h>
0026 #include <trackbase/TrkrDefs.h>
0027 #include <trackbase/TrkrHit.h>
0028 #include <trackbase/TrkrHitSet.h>
0029 #include <trackbase/TrkrHitSetContainer.h>
0030
0031 #include <phgarfield/PHGarfield.h>
0032 #include <TPolyLine3D.h>
0033
0034 #include <g4detectors/PHG4TpcGeom.h>
0035 #include <g4detectors/PHG4TpcGeomContainer.h>
0036
0037 #include <algorithm>
0038 #include <array>
0039 #include <cmath>
0040 #include <iostream>
0041 #include <limits>
0042 #include <map>
0043 #include <set>
0044 #include <utility>
0045 #include <vector>
0046
0047 namespace
0048 {
0049 constexpr unsigned int FirstLayer = 7;
0050 constexpr unsigned int LastLayer = 54;
0051 constexpr unsigned int NLayers = LastLayer - FirstLayer + 1;
0052 constexpr unsigned int NSides = 2;
0053 constexpr unsigned int NSectors = 12;
0054 constexpr double PhiConsistencyTolerance = 1.0e-10;
0055 constexpr unsigned char InvalidConfidenceTier = std::numeric_limits<unsigned char>::max();
0056
0057 double wrap_phi(double phi)
0058 {
0059 while (phi > M_PI) { phi -= 2.0 * M_PI;
0060 }
0061 while (phi <= -M_PI) { phi += 2.0 * M_PI;
0062 }
0063 return phi;
0064 }
0065
0066 double unwrap_phi_near(double phi, const double reference)
0067 {
0068 while (phi - reference > M_PI) { phi -= 2.0 * M_PI;
0069 }
0070 while (phi - reference <= -M_PI) { phi += 2.0 * M_PI;
0071 }
0072 return phi;
0073 }
0074
0075 double clamp_unit(const double value)
0076 {
0077 return std::max(0.0, std::min(1.0, value));
0078 }
0079
0080 double phi_sample_fraction(const unsigned int sample)
0081 {
0082 if (TpcCrossingFinder::NPhiSamples <= 1U) { return 0.0;
0083 }
0084 return static_cast<double>(sample) / static_cast<double>(TpcCrossingFinder::NPhiSamples - 1U);
0085 }
0086
0087 float finite_float_or_nan(const double value)
0088 {
0089 return std::isfinite(value) ? static_cast<float>(value) : std::numeric_limits<float>::quiet_NaN();
0090 }
0091 }
0092
0093 TpcCrossingFinder::TpcCrossingFinder(const std::string& name)
0094 : SubsysReco(name)
0095 {
0096 }
0097
0098 TpcCrossingFinder::~TpcCrossingFinder()
0099 {
0100 delete m_idealPadMap;
0101 m_idealPadMap = nullptr;
0102 delete m_garfield;
0103 m_garfield = nullptr;
0104 }
0105
0106 int TpcCrossingFinder::InitRun(PHCompositeNode* topNode)
0107 {
0108 if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK) { return Fun4AllReturnCodes::ABORTRUN;
0109 }
0110 if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK) { return Fun4AllReturnCodes::ABORTRUN;
0111 }
0112
0113 delete m_idealPadMap;
0114 m_idealPadMap = new IdealPadMap();
0115 if (m_idealPadMap->load_from_cdb(Verbosity()) != 0 || !m_idealPadMap->is_loaded())
0116 {
0117 std::cerr << Name() << "::InitRun - failed to load IdealPadMap" << std::endl;
0118 return Fun4AllReturnCodes::ABORTRUN;
0119 }
0120
0121 PHG4TpcGeom* layergeom = m_geomContainerTpc->GetLayerCellGeom(20);
0122 if (layergeom)
0123 {
0124 const double rot_x = layergeom->get_rot_x();
0125 const double rot_y = layergeom->get_rot_y();
0126 const double rot_z = layergeom->get_rot_z();
0127 const double place_x = layergeom->get_place_x();
0128 const double place_y = layergeom->get_place_y();
0129 const double place_z = layergeom->get_place_z();
0130 if (use_survey_geometry)
0131 {
0132 m_tpcMove = {{place_x, place_y, place_z}};
0133 m_tpcRotations = {{{{rot_x, rot_y, rot_z}}, {{0.0, 0.0, 0.0}}}};
0134 }
0135 }
0136
0137 delete m_garfield;
0138 const std::string electricFieldMap = CDBInterface::instance()->getUrl("Tpc_PolySeeding_EField");
0139
0140 const auto kefffile = CDBInterface::instance()->getUrl("Tpc_PolyClusterizer_kEff");
0141
0142 if (!kefffile.empty())
0143 {
0144 auto keffcdbtree = std::make_unique<CDBTTree>(kefffile);
0145 keffcdbtree->LoadCalibrations();
0146 m_kEffSide0 = keffcdbtree->GetSingleFloatValue("keffside0");
0147 m_kEffSide1 = keffcdbtree->GetSingleFloatValue("keffside1");
0148 }
0149
0150 m_garfield = new PHGarfield(Name() + "_PHGarfield", electricFieldMap, m_kEffSide0, m_kEffSide1);
0151 configure_garfield(m_garfield);
0152 if (m_garfield->InitRun(topNode) != Fun4AllReturnCodes::EVENT_OK)
0153 {
0154 std::cerr << Name() << "::InitRun - PHGarfield InitRun failed" << std::endl;
0155 return Fun4AllReturnCodes::ABORTRUN;
0156 }
0157
0158 if (!build_drift_lookup()) { return Fun4AllReturnCodes::ABORTRUN;
0159 }
0160
0161 m_event = 0;
0162 return Fun4AllReturnCodes::EVENT_OK;
0163 }
0164
0165 int TpcCrossingFinder::getNodes(PHCompositeNode* topNode)
0166 {
0167 m_assembledTracks = findNode::getClass<Tpc_AssembledTrackContainer>(topNode, m_inputNodeName);
0168 if (!m_assembledTracks)
0169 {
0170 std::cerr << Name() << "::getNodes - missing " << m_inputNodeName << std::endl;
0171 return Fun4AllReturnCodes::ABORTRUN;
0172 }
0173
0174 m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0175 if (!m_hits)
0176 {
0177 std::cerr << Name() << "::getNodes - missing TRKR_HITSET" << std::endl;
0178 return Fun4AllReturnCodes::ABORTRUN;
0179 }
0180
0181 m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0182 if (!m_clusterMap && Verbosity() > 0)
0183 {
0184 std::cout << Name() << "::getNodes - optional TRKR_CLUSTER node not found" << std::endl;
0185 }
0186
0187 m_geomContainerTpc = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0188 if (!m_geomContainerTpc)
0189 {
0190 std::cerr << Name() << "::getNodes - missing TPCGEOMCONTAINER" << std::endl;
0191 return Fun4AllReturnCodes::ABORTRUN;
0192 }
0193
0194 m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, m_vertexMapNodeName);
0195 if (!m_vertexMap && Verbosity() > 0)
0196 {
0197 std::cout << Name() << "::getNodes - optional " << m_vertexMapNodeName << " node not found" << std::endl;
0198 }
0199
0200 return Fun4AllReturnCodes::EVENT_OK;
0201 }
0202
0203 int TpcCrossingFinder::createNodes(PHCompositeNode* topNode)
0204 {
0205 PHNodeIterator iter(topNode);
0206 PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0207 if (!dstNode)
0208 {
0209 dstNode = new PHCompositeNode("DST");
0210 topNode->addNode(dstNode);
0211 }
0212
0213 m_decisions = findNode::getClass<TpcCrossingDecisionContainerv1>(topNode, m_outputNodeName);
0214 if (!m_decisions)
0215 {
0216 m_decisions = new TpcCrossingDecisionContainerv1();
0217 PHIODataNode<PHObject>* node = new PHIODataNode<PHObject>(m_decisions, m_outputNodeName, "PHObject");
0218 dstNode->addNode(node);
0219 std::cout << Name() << "::createNodes - created " << m_outputNodeName << " node" << std::endl;
0220 }
0221
0222 return Fun4AllReturnCodes::EVENT_OK;
0223 }
0224
0225 void TpcCrossingFinder::configure_garfield(PHGarfield* garfield) const
0226 {
0227 if (!garfield) { return;
0228 }
0229
0230 garfield->MoveTpc(m_tpcMove[0], m_tpcMove[1], m_tpcMove[2]);
0231 for (const auto& rotation : m_tpcRotations)
0232 {
0233 garfield->RotateTpc(rotation[0], rotation[1], rotation[2]);
0234 }
0235 garfield->SetCMVoltageDefault(m_cmVoltageDefault);
0236 }
0237
0238 unsigned int TpcCrossingFinder::drift_lookup_index(const unsigned int layer_index,
0239 const unsigned int side,
0240 const unsigned int sector,
0241 const unsigned int sample)
0242 {
0243 return (((layer_index * NSides + side) * NSectors + sector) * NPhiSamples + sample);
0244 }
0245
0246 bool TpcCrossingFinder::build_drift_lookup()
0247 {
0248 if (!m_idealPadMap || !m_garfield) { return false;
0249 }
0250 if (m_reverseDriftStepNs <= 0.0 || !std::isfinite(m_reverseDriftStepNs)) { return false;
0251 }
0252
0253 m_maxLookupTimeNs = std::numeric_limits<double>::max();
0254 for (DriftPolyline& polyline : m_driftLookup)
0255 {
0256 polyline.phi = 0.0;
0257 polyline.points.clear();
0258 }
0259
0260 unsigned int nbuilt = 0;
0261 for (unsigned int layer = FirstLayer; layer <= LastLayer; ++layer)
0262 {
0263 const unsigned int layer_index = layer - FirstLayer;
0264 const double radius = m_idealPadMap->get_radius(layer);
0265 const unsigned int pads_per_sector = m_idealPadMap->get_pads_per_sector_for_layer(layer);
0266 if (!std::isfinite(radius) || pads_per_sector == 0U) { return false;
0267 }
0268
0269 for (unsigned int side = 0; side < NSides; ++side)
0270 {
0271 const double z0 = (side == 0U) ? m_startZSouth : m_startZNorth;
0272 for (unsigned int sector = 0; sector < NSectors; ++sector)
0273 {
0274 for (unsigned int sample = 0; sample < NPhiSamples; ++sample)
0275 {
0276 const unsigned int local_phibin = static_cast<unsigned int>(std::llround(
0277 phi_sample_fraction(sample) * static_cast<double>(pads_per_sector - 1U)));
0278 const unsigned int global_pad = sector * pads_per_sector + local_phibin;
0279 const double phi_local = m_idealPadMap->get_phi(side, sector, layer, local_phibin);
0280 const double phi_global = m_idealPadMap->get_phi(side, layer, global_pad);
0281 if (!std::isfinite(phi_local) || !std::isfinite(phi_global)) { return false;
0282 }
0283 if (std::abs(wrap_phi(phi_global - phi_local)) > PhiConsistencyTolerance) { return false;
0284 }
0285
0286 const double x0 = radius * std::cos(phi_local);
0287 const double y0 = radius * std::sin(phi_local);
0288 TPolyLine3D* drift = m_garfield->ReverseDrift(x0, y0, z0, m_reverseDriftStepNs);
0289 if (!drift || drift->GetN() <= 0)
0290 {
0291 delete drift;
0292 return false;
0293 }
0294
0295 const int npoints = drift->GetN();
0296 const Float_t* xyz = drift->GetP();
0297 if (!xyz || npoints <= 0)
0298 {
0299 delete drift;
0300 return false;
0301 }
0302
0303 DriftPolyline& polyline = m_driftLookup[drift_lookup_index(layer_index, side, sector, sample)];
0304 polyline.phi = phi_local;
0305 polyline.points.resize(static_cast<std::size_t>(npoints));
0306 for (int ipoint = 0; ipoint < npoints; ++ipoint)
0307 {
0308 const int idx = 3 * ipoint;
0309 DriftPoint& point = polyline.points[static_cast<std::size_t>(ipoint)];
0310 const double output_r = std::hypot(static_cast<double>(xyz[idx]), static_cast<double>(xyz[idx + 1]));
0311 const double output_phi = unwrap_phi_near(std::atan2(static_cast<double>(xyz[idx + 1]), static_cast<double>(xyz[idx])), phi_local);
0312 point.delta_r = static_cast<float>(output_r - radius);
0313 point.delta_phi = static_cast<float>(output_phi - phi_local);
0314 point.z = xyz[idx + 2];
0315 }
0316 m_maxLookupTimeNs = std::min(m_maxLookupTimeNs, static_cast<double>(npoints - 1) * m_reverseDriftStepNs);
0317 delete drift;
0318 ++nbuilt;
0319 }
0320 }
0321 }
0322 }
0323
0324 if (Verbosity() > 0)
0325 {
0326 std::cout << Name() << "::build_drift_lookup - built " << nbuilt
0327 << " drift polylines max_lookup_time_ns=" << m_maxLookupTimeNs << std::endl;
0328 }
0329 return nbuilt == NLayers * NSides * NSectors * NPhiSamples && std::isfinite(m_maxLookupTimeNs) && m_maxLookupTimeNs > 0.0;
0330 }
0331
0332 bool TpcCrossingFinder::sample_drift_lookup(const unsigned int layer,
0333 const unsigned int side,
0334 const unsigned int pad,
0335 const unsigned int tbin,
0336 const short crossing,
0337 double& x,
0338 double& y,
0339 double& z) const
0340 {
0341 if (!m_idealPadMap) { return false;
0342 }
0343 if (layer < FirstLayer || layer > LastLayer) { return false;
0344 }
0345 if (side >= NSides) { return false;
0346 }
0347 if (m_reverseDriftStepNs <= 0.0 || !std::isfinite(m_reverseDriftStepNs)) { return false;
0348 }
0349
0350 const unsigned int pads_per_sector = m_idealPadMap->get_pads_per_sector_for_layer(layer);
0351 if (pads_per_sector == 0U) { return false;
0352 }
0353 const unsigned int sector = pad / pads_per_sector;
0354 if (sector >= NSectors) { return false;
0355 }
0356
0357 const double hit_radius = m_idealPadMap->get_radius(layer);
0358 const double hit_phi = m_idealPadMap->get_phi(side, layer, pad);
0359 if (!std::isfinite(hit_radius) || !std::isfinite(hit_phi)) { return false;
0360 }
0361
0362 const double target_time_ns = (static_cast<double>(tbin) - m_t0) * m_tpcAdcClock
0363 - static_cast<double>(crossing) * m_crossingPeriodNs;
0364 if (target_time_ns <= 0.0 || !std::isfinite(target_time_ns)) { return false;
0365 }
0366
0367 const unsigned int layer_index = layer - FirstLayer;
0368 std::array<const DriftPolyline*, NPhiSamples> samples{};
0369 std::array<double, NPhiSamples> sample_phi{};
0370 for (unsigned int sample = 0; sample < NPhiSamples; ++sample)
0371 {
0372 samples[sample] = &m_driftLookup[drift_lookup_index(layer_index, side, sector, sample)];
0373 if (samples[sample]->points.empty()) { return false;
0374 }
0375 sample_phi[sample] = samples[sample]->phi;
0376 if (sample > 0U) { sample_phi[sample] = unwrap_phi_near(sample_phi[sample], sample_phi[sample - 1U]);
0377 }
0378 }
0379
0380 const double unwrapped_hit_phi = unwrap_phi_near(hit_phi, sample_phi[NPhiSamples / 2U]);
0381 const bool increasing = sample_phi[NPhiSamples - 1U] >= sample_phi[0];
0382
0383 bool bracket_found = false;
0384 unsigned int sample0 = 0;
0385 unsigned int sample1 = 0;
0386 double phi_fraction = 0.0;
0387 if ((increasing && unwrapped_hit_phi <= sample_phi[0]) || (!increasing && unwrapped_hit_phi >= sample_phi[0]))
0388 {
0389 bracket_found = true;
0390 }
0391 else if ((increasing && unwrapped_hit_phi >= sample_phi[NPhiSamples - 1U]) ||
0392 (!increasing && unwrapped_hit_phi <= sample_phi[NPhiSamples - 1U]))
0393 {
0394 sample0 = NPhiSamples - 1U;
0395 sample1 = NPhiSamples - 1U;
0396 bracket_found = true;
0397 }
0398 else
0399 {
0400 for (unsigned int sample = 0; sample + 1U < NPhiSamples; ++sample)
0401 {
0402 const bool in_interval = increasing ?
0403 (unwrapped_hit_phi >= sample_phi[sample] && unwrapped_hit_phi <= sample_phi[sample + 1U]) :
0404 (unwrapped_hit_phi <= sample_phi[sample] && unwrapped_hit_phi >= sample_phi[sample + 1U]);
0405 if (!in_interval) { continue;
0406 }
0407
0408 sample0 = sample;
0409 sample1 = sample + 1U;
0410 const double denom = sample_phi[sample1] - sample_phi[sample0];
0411 phi_fraction = (denom != 0.0) ? clamp_unit((unwrapped_hit_phi - sample_phi[sample0]) / denom) : 0.0;
0412 bracket_found = true;
0413 break;
0414 }
0415 }
0416 if (!bracket_found) { return false;
0417 }
0418
0419 auto sample_time = [this, target_time_ns](const DriftPolyline& polyline,
0420 double& delta_r,
0421 double& delta_phi,
0422 double& point_z) -> bool
0423 {
0424 const int npoints = static_cast<int>(polyline.points.size());
0425 if (npoints <= 0) { return false;
0426 }
0427 const double max_time_ns = static_cast<double>(npoints - 1) * m_reverseDriftStepNs;
0428 if (target_time_ns > max_time_ns) { return false;
0429 }
0430
0431 const double fbin = target_time_ns / m_reverseDriftStepNs;
0432 const int i0 = std::min(static_cast<int>(std::floor(fbin)), npoints - 1);
0433 const int i1 = std::min(i0 + 1, npoints - 1);
0434 const double frac = fbin - static_cast<double>(i0);
0435
0436 const DriftPoint& p0 = polyline.points[static_cast<std::size_t>(i0)];
0437 const DriftPoint& p1 = polyline.points[static_cast<std::size_t>(i1)];
0438 const double dphi0 = static_cast<double>(p0.delta_phi);
0439 const double dphi1 = unwrap_phi_near(static_cast<double>(p1.delta_phi), dphi0);
0440 delta_r = static_cast<double>(p0.delta_r) + frac * static_cast<double>(p1.delta_r - p0.delta_r);
0441 delta_phi = dphi0 + frac * (dphi1 - dphi0);
0442 point_z = static_cast<double>(p0.z) + frac * static_cast<double>(p1.z - p0.z);
0443 return std::isfinite(delta_r) && std::isfinite(delta_phi) && std::isfinite(point_z);
0444 };
0445
0446 double delta_r0 = 0.0;
0447 double delta_phi0 = 0.0;
0448 double z0 = 0.0;
0449 const bool valid0 = sample_time(*samples[sample0], delta_r0, delta_phi0, z0);
0450 double delta_r1 = 0.0;
0451 double delta_phi1 = 0.0;
0452 double z1 = 0.0;
0453 const bool same_sample = sample0 == sample1;
0454 const bool valid1 = same_sample ? valid0 : sample_time(*samples[sample1], delta_r1, delta_phi1, z1);
0455 if (same_sample)
0456 {
0457 delta_r1 = delta_r0;
0458 delta_phi1 = delta_phi0;
0459 z1 = z0;
0460 }
0461 if (!valid0 && !valid1) { return false;
0462 }
0463
0464 double delta_r = 0.0;
0465 double delta_phi = 0.0;
0466 if (valid0 && valid1 && !same_sample)
0467 {
0468 delta_r = delta_r0 + phi_fraction * (delta_r1 - delta_r0);
0469 const double unwrapped_delta_phi1 = unwrap_phi_near(delta_phi1, delta_phi0);
0470 delta_phi = delta_phi0 + phi_fraction * (unwrapped_delta_phi1 - delta_phi0);
0471 z = z0 + phi_fraction * (z1 - z0);
0472 }
0473 else if (valid0)
0474 {
0475 delta_r = delta_r0;
0476 delta_phi = delta_phi0;
0477 z = z0;
0478 }
0479 else
0480 {
0481 delta_r = delta_r1;
0482 delta_phi = delta_phi1;
0483 z = z1;
0484 }
0485
0486 const double radius = hit_radius + delta_r;
0487 const double output_phi = unwrapped_hit_phi + delta_phi;
0488 x = radius * std::cos(output_phi);
0489 y = radius * std::sin(output_phi);
0490 return std::isfinite(x) && std::isfinite(y) && std::isfinite(z);
0491 }
0492
0493 bool TpcCrossingFinder::make_xyz_point(TrkrDefs::hitsetkey hsk,
0494 TrkrDefs::hitkey hk,
0495 short crossing,
0496 Point& p) const
0497 {
0498 const unsigned int layer = TrkrDefs::getLayer(hsk);
0499 const unsigned int hit_side = TpcDefs::getSide(hsk);
0500 const unsigned int pad = TpcDefs::getPad(hk);
0501 const unsigned int tbin = TpcDefs::getTBin(hk);
0502
0503 double x = 0.0;
0504 double y = 0.0;
0505 double z = 0.0;
0506 if (!sample_drift_lookup(layer, hit_side, pad, tbin, crossing, x, y, z)) { return false;
0507 }
0508
0509 p.hitsetkey = hsk;
0510 p.hitkey = hk;
0511 p.layer = layer;
0512 p.side = hit_side;
0513 p.pad = pad;
0514 p.tbin = tbin;
0515 p.x = x;
0516 p.y = y;
0517 p.z = z;
0518 return true;
0519 }
0520
0521 bool TpcCrossingFinder::find_time_extrema(const Tpc_AssembledTrack* track,
0522 TrkrDefs::hitsetkey& min_hsk,
0523 TrkrDefs::hitkey& min_hk,
0524 TrkrDefs::hitsetkey& max_hsk,
0525 TrkrDefs::hitkey& max_hk) const
0526 {
0527 if (!track || track->size_hit_indices() == 0) { return false;
0528 }
0529
0530 unsigned int min_tbin = std::numeric_limits<unsigned int>::max();
0531 unsigned int max_tbin = 0;
0532 bool found = false;
0533 for (unsigned int ih = 0; ih < track->size_hit_indices(); ++ih)
0534 {
0535 const Tpc_AssembledTrack::HitIndex hi = track->get_hit_index(ih);
0536 const unsigned int tbin = TpcDefs::getTBin(hi.second);
0537 if (!found || tbin < min_tbin)
0538 {
0539 min_tbin = tbin;
0540 min_hsk = hi.first;
0541 min_hk = hi.second;
0542 }
0543 if (!found || tbin > max_tbin)
0544 {
0545 max_tbin = tbin;
0546 max_hsk = hi.first;
0547 max_hk = hi.second;
0548 }
0549 found = true;
0550 }
0551 return found;
0552 }
0553
0554 std::set<short> TpcCrossingFinder::get_available_crossings() const
0555 {
0556 std::set<short> available_crossings;
0557 if (!m_vertexMap) { return available_crossings;
0558 }
0559
0560 for (auto & iter : *m_vertexMap)
0561 {
0562 const SvtxVertex* vertex = iter.second;
0563 if (!vertex) { continue;
0564 }
0565 available_crossings.insert(vertex->get_beam_crossing());
0566 }
0567 return available_crossings;
0568 }
0569
0570 std::set<short> TpcCrossingFinder::get_intt_crossings() const
0571 {
0572 std::set<short> intt_crossings;
0573 if (!m_clusterMap) { return intt_crossings;
0574 }
0575
0576 for (const auto& hitsetkey : m_clusterMap->getHitSetKeys())
0577 {
0578 if (static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(hitsetkey)) != TrkrDefs::TrkrId::inttId) { continue;
0579 }
0580 intt_crossings.insert(static_cast<short>(InttDefs::getTimeBucketId(hitsetkey)));
0581 }
0582 return intt_crossings;
0583 }
0584
0585 std::map<short, std::vector<TpcCrossingFinder::SiliconVertexHypothesis>>
0586 TpcCrossingFinder::get_vertices_by_crossing() const
0587 {
0588 std::map<short, std::vector<SiliconVertexHypothesis>> vertices_by_crossing;
0589 if (!m_vertexMap) { return vertices_by_crossing;
0590 }
0591
0592 for (auto & iter : *m_vertexMap)
0593 {
0594 const SvtxVertex* vertex = iter.second;
0595 if (!vertex) { continue;
0596 }
0597
0598 SiliconVertexHypothesis hyp;
0599 hyp.crossing = vertex->get_beam_crossing();
0600 hyp.vertex_id = iter.first;
0601 hyp.x = vertex->get_x();
0602 hyp.y = vertex->get_y();
0603 hyp.z = vertex->get_z();
0604 const double errzz = vertex->get_error(2, 2);
0605 hyp.sigma_z = errzz >= 0.0 ? std::sqrt(errzz) : std::numeric_limits<double>::quiet_NaN();
0606 hyp.ntracks = vertex->size_tracks();
0607 vertices_by_crossing[hyp.crossing].push_back(hyp);
0608 }
0609 return vertices_by_crossing;
0610 }
0611
0612 std::vector<std::pair<TrkrDefs::hitsetkey, TrkrDefs::hitkey>>
0613 TpcCrossingFinder::select_representatives(const Tpc_AssembledTrack* track,
0614 TrkrDefs::hitsetkey min_hsk,
0615 TrkrDefs::hitkey min_hk,
0616 TrkrDefs::hitsetkey max_hsk,
0617 TrkrDefs::hitkey max_hk) const
0618 {
0619 std::vector<std::pair<TrkrDefs::hitsetkey, TrkrDefs::hitkey>> representatives;
0620 auto add_unique = [&representatives](const Tpc_AssembledTrack::HitIndex& hi)
0621 {
0622 if (std::find(representatives.begin(), representatives.end(), hi) == representatives.end())
0623 {
0624 representatives.push_back(hi);
0625 }
0626 };
0627
0628 add_unique({min_hsk, min_hk});
0629 add_unique({max_hsk, max_hk});
0630 if (!track) { return representatives;
0631 }
0632
0633 std::map<unsigned int, Tpc_AssembledTrack::HitIndex> hit_by_layer;
0634 for (unsigned int ih = 0; ih < track->size_hit_indices(); ++ih)
0635 {
0636 const Tpc_AssembledTrack::HitIndex hi = track->get_hit_index(ih);
0637 hit_by_layer.emplace(TrkrDefs::getLayer(hi.first), hi);
0638 }
0639 if (hit_by_layer.empty()) { return representatives;
0640 }
0641
0642 add_unique(hit_by_layer.begin()->second);
0643 add_unique(hit_by_layer.rbegin()->second);
0644
0645 std::vector<unsigned int> layers;
0646 layers.reserve(hit_by_layer.size());
0647 for (const auto& item : hit_by_layer) { layers.push_back(item.first);
0648 }
0649 const unsigned int fractions[] = {1U, 2U};
0650 for (const unsigned int fraction : fractions)
0651 {
0652 if (layers.size() < 3U) { continue;
0653 }
0654 const std::size_t index = (fraction * (layers.size() - 1U)) / 3U;
0655 add_unique(hit_by_layer[layers[index]]);
0656 }
0657
0658 return representatives;
0659 }
0660
0661 TpcCrossingFinder::ZFitResult TpcCrossingFinder::estimate_tpc_z0_diagnostics(std::vector<Point> points) const
0662 {
0663 ZFitResult result;
0664 if (points.size() < 2U) { return result;
0665 }
0666 std::sort(points.begin(), points.end(), [](const Point& a, const Point& b) {
0667 return a.layer < b.layer;
0668 });
0669
0670 std::vector<double> s(points.size(), 0.0);
0671 for (std::size_t i = 1; i < points.size(); ++i)
0672 {
0673 s[i] = s[i - 1] + std::hypot(points[i].x - points[i - 1].x, points[i].y - points[i - 1].y);
0674 }
0675
0676 std::vector<Tpc_FittingTools::FitPoint> fit_points;
0677 fit_points.reserve(points.size());
0678 std::vector<Tpc_FittingTools::FitPoint> radial_fit_points;
0679 radial_fit_points.reserve(points.size());
0680 for (std::size_t i = 0; i < points.size(); ++i)
0681 {
0682 fit_points.emplace_back(s[i], points[i].z, 1.0);
0683 radial_fit_points.emplace_back(std::hypot(points[i].x, points[i].y), points[i].z, 1.0);
0684 }
0685 const Tpc_FittingTools::LineFit zfit = Tpc_FittingTools::fitLine(fit_points);
0686 const Tpc_FittingTools::LineFit radial_zfit = Tpc_FittingTools::fitLine(radial_fit_points);
0687 if (!zfit.ok || !radial_zfit.ok) { return result;
0688 }
0689
0690 double best_s = s.front();
0691 double best_d2 = std::numeric_limits<double>::max();
0692 for (std::size_t i = 0; i + 1U < points.size(); ++i)
0693 {
0694 const double vx = points[i + 1U].x - points[i].x;
0695 const double vy = points[i + 1U].y - points[i].y;
0696 const double len2 = vx * vx + vy * vy;
0697 double frac = 0.0;
0698 if (len2 > 0.0)
0699 {
0700 frac = clamp_unit(-(points[i].x * vx + points[i].y * vy) / len2);
0701 }
0702 const double x = points[i].x + frac * vx;
0703 const double y = points[i].y + frac * vy;
0704 const double d2 = x * x + y * y;
0705 if (d2 < best_d2)
0706 {
0707 best_d2 = d2;
0708 best_s = s[i] + frac * (s[i + 1U] - s[i]);
0709 }
0710 }
0711
0712 result.valid = std::isfinite(best_s) && std::isfinite(best_d2);
0713 result.slope = zfit.slope;
0714 result.intercept = zfit.intercept;
0715 result.chi2 = zfit.chi2;
0716 result.ndf = zfit.ndof;
0717 result.s_at_pca = best_s;
0718 result.minimum_radius = std::sqrt(best_d2);
0719 result.z_at_pca = zfit.slope * best_s + zfit.intercept;
0720 result.z_at_r0 = radial_zfit.intercept;
0721 result.points = std::move(points);
0722 result.path_length.reserve(s.size());
0723 for (const double value : s) { result.path_length.push_back(finite_float_or_nan(value));
0724 }
0725 result.valid = result.valid && std::isfinite(result.z_at_pca) && std::isfinite(result.z_at_r0);
0726 return result;
0727 }
0728
0729 bool TpcCrossingFinder::estimate_tpc_z0(std::vector<Point>& points, double& z0) const
0730 {
0731 const ZFitResult result = estimate_tpc_z0_diagnostics(points);
0732 if (!result.valid) { return false;
0733 }
0734 z0 = result.z_at_r0;
0735 return true;
0736 }
0737
0738 bool TpcCrossingFinder::point_in_tpc(const Point& p) const
0739 {
0740 const double r = std::hypot(p.x, p.y);
0741 const double inner = m_idealPadMap ? m_idealPadMap->get_radius(FirstLayer) : 0.0;
0742 const double outer = m_idealPadMap ? m_idealPadMap->get_radius(LastLayer) : 0.0;
0743 return std::isfinite(r) && std::isfinite(p.z) &&
0744 r >= inner - m_radialTolerance &&
0745 r <= outer + m_radialTolerance &&
0746 p.z >= -m_tpcHalfLength - m_zTolerance &&
0747 p.z <= m_tpcHalfLength + m_zTolerance;
0748 }
0749
0750 bool TpcCrossingFinder::point_in_correct_side(const Point& p) const
0751 {
0752 if (p.side == 0U) { return p.z >= -m_tpcHalfLength - m_zTolerance && p.z <= m_centralMembraneTolerance;
0753 }
0754 if (p.side == 1U) { return p.z >= -m_centralMembraneTolerance && p.z <= m_tpcHalfLength + m_zTolerance;
0755 }
0756 return false;
0757 }
0758
0759 TpcCrossingFinder::Candidate
0760 TpcCrossingFinder::test_candidate(const Tpc_AssembledTrack* track,
0761 short crossing,
0762 TrkrDefs::hitsetkey min_hsk,
0763 TrkrDefs::hitkey min_hk,
0764 TrkrDefs::hitsetkey max_hsk,
0765 TrkrDefs::hitkey max_hk,
0766 const std::map<short, std::vector<SiliconVertexHypothesis>>& vertices_by_crossing) const
0767 {
0768 Candidate candidate;
0769 candidate.crossing = crossing;
0770 candidate.silicon_vertex_id = std::numeric_limits<unsigned int>::max();
0771 candidate.tpc_z0 = std::numeric_limits<double>::quiet_NaN();
0772 candidate.silicon_vertex_z = std::numeric_limits<double>::quiet_NaN();
0773 candidate.delta_z = std::numeric_limits<double>::quiet_NaN();
0774
0775 TpcCrossingCandidate& qa = candidate.qa;
0776 qa.crossing = crossing;
0777 qa.is_available_from_intt = true;
0778 qa.passes_time_window = true;
0779 qa.was_tested = true;
0780 qa.rejection_status = static_cast<unsigned char>(TpcCrossingStatus::Unknown);
0781 qa.max_lookup_time_ns = finite_float_or_nan(m_maxLookupTimeNs);
0782 qa.min_tbin_time_crossing0_ns = finite_float_or_nan((static_cast<double>(TpcDefs::getTBin(min_hk)) - m_t0) * m_tpcAdcClock);
0783 qa.max_tbin_time_crossing0_ns = finite_float_or_nan((static_cast<double>(TpcDefs::getTBin(max_hk)) - m_t0) * m_tpcAdcClock);
0784 qa.candidate_min_time_ns = finite_float_or_nan(static_cast<double>(qa.min_tbin_time_crossing0_ns) - static_cast<double>(crossing) * m_crossingPeriodNs);
0785 qa.candidate_max_time_ns = finite_float_or_nan(static_cast<double>(qa.max_tbin_time_crossing0_ns) - static_cast<double>(crossing) * m_crossingPeriodNs);
0786 qa.min_time_margin_ns = qa.candidate_min_time_ns;
0787 qa.max_time_margin_ns = finite_float_or_nan(m_maxLookupTimeNs - static_cast<double>(qa.candidate_max_time_ns));
0788 qa.candidate_qa_bits |= FromInttCrossing | PassedTimeWindow;
0789
0790 auto fill_point_summary = [](const Point& point,
0791 unsigned int& layer,
0792 unsigned int& side,
0793 unsigned int& pad,
0794 unsigned int& tbin,
0795 float& x,
0796 float& y,
0797 float& z,
0798 float& r,
0799 float& phi)
0800 {
0801 layer = point.layer;
0802 side = point.side;
0803 pad = point.pad;
0804 tbin = point.tbin;
0805 x = finite_float_or_nan(point.x);
0806 y = finite_float_or_nan(point.y);
0807 z = finite_float_or_nan(point.z);
0808 r = finite_float_or_nan(std::hypot(point.x, point.y));
0809 phi = finite_float_or_nan(std::atan2(point.y, point.x));
0810 };
0811
0812 auto distance_to_padplane = [this](const Point& point) -> float
0813 {
0814 const double padplane_z = point.side == 0U ? -m_tpcHalfLength : m_tpcHalfLength;
0815 return finite_float_or_nan(std::abs(point.z - padplane_z));
0816 };
0817
0818 Point min_point;
0819 Point max_point;
0820 qa.min_endpoint_garfield_valid = make_xyz_point(min_hsk, min_hk, crossing, min_point);
0821 qa.max_endpoint_garfield_valid = make_xyz_point(max_hsk, max_hk, crossing, max_point);
0822 if (qa.min_endpoint_garfield_valid)
0823 {
0824 fill_point_summary(min_point, qa.min_time_layer, qa.min_time_side, qa.min_time_pad, qa.min_time_tbin,
0825 qa.min_time_x, qa.min_time_y, qa.min_time_z, qa.min_time_r, qa.min_time_phi);
0826 qa.min_time_distance_to_central_membrane = finite_float_or_nan(std::abs(min_point.z));
0827 qa.min_time_distance_to_padplane = distance_to_padplane(min_point);
0828 qa.min_time_inside_tpc = point_in_tpc(min_point);
0829 qa.min_time_correct_side = point_in_correct_side(min_point);
0830 qa.candidate_qa_bits |= MinEndpointGarfieldOK;
0831 }
0832 if (qa.max_endpoint_garfield_valid)
0833 {
0834 fill_point_summary(max_point, qa.max_time_layer, qa.max_time_side, qa.max_time_pad, qa.max_time_tbin,
0835 qa.max_time_x, qa.max_time_y, qa.max_time_z, qa.max_time_r, qa.max_time_phi);
0836 qa.max_time_distance_to_central_membrane = finite_float_or_nan(std::abs(max_point.z));
0837 qa.max_time_distance_to_padplane = distance_to_padplane(max_point);
0838 qa.max_time_inside_tpc = point_in_tpc(max_point);
0839 qa.max_time_correct_side = point_in_correct_side(max_point);
0840 qa.candidate_qa_bits |= MaxEndpointGarfieldOK;
0841 }
0842
0843 if (!qa.min_endpoint_garfield_valid || !qa.max_endpoint_garfield_valid)
0844 {
0845 qa.first_failed_stage = static_cast<unsigned char>(TpcCrossingCandidateStage::GarfieldValid);
0846 candidate.rejection_status = static_cast<unsigned char>(TpcCrossingStatus::GarfieldInvalid);
0847 qa.rejection_status = candidate.rejection_status;
0848 return candidate;
0849 }
0850
0851 qa.endpoints_inside_tpc = qa.min_time_inside_tpc && qa.max_time_inside_tpc;
0852 if (qa.endpoints_inside_tpc) { qa.candidate_qa_bits |= EndpointsInsideTPC;
0853 }
0854 if (!qa.endpoints_inside_tpc)
0855 {
0856 qa.first_failed_stage = static_cast<unsigned char>(TpcCrossingCandidateStage::InsideTpc);
0857 candidate.rejection_status = static_cast<unsigned char>(TpcCrossingStatus::OutsideTpcVolume);
0858 qa.rejection_status = candidate.rejection_status;
0859 return candidate;
0860 }
0861
0862 qa.endpoints_on_correct_side = qa.min_time_correct_side && qa.max_time_correct_side;
0863 if (qa.endpoints_on_correct_side) { qa.candidate_qa_bits |= EndpointsCorrectSide;
0864 }
0865 if (!qa.endpoints_on_correct_side)
0866 {
0867 qa.first_failed_stage = static_cast<unsigned char>(TpcCrossingCandidateStage::CorrectSide);
0868 candidate.rejection_status = static_cast<unsigned char>(TpcCrossingStatus::WrongTpcSide);
0869 qa.rejection_status = candidate.rejection_status;
0870 return candidate;
0871 }
0872
0873 std::vector<Point> points;
0874 for (const auto& hi : select_representatives(track, min_hsk, min_hk, max_hsk, max_hk))
0875 {
0876 Point p;
0877 if (!make_xyz_point(hi.first, hi.second, crossing, p)) { continue;
0878 }
0879 if (!point_in_tpc(p) || !point_in_correct_side(p)) { continue;
0880 }
0881 points.push_back(p);
0882 }
0883
0884 const ZFitResult fit = estimate_tpc_z0_diagnostics(points);
0885 qa.fit_valid = fit.valid;
0886 if (!fit.valid)
0887 {
0888 qa.first_failed_stage = static_cast<unsigned char>(TpcCrossingCandidateStage::FitValid);
0889 candidate.rejection_status = static_cast<unsigned char>(TpcCrossingStatus::FitFailed);
0890 qa.rejection_status = candidate.rejection_status;
0891 return candidate;
0892 }
0893
0894 qa.candidate_qa_bits |= FitSuccessful;
0895 qa.n_fit_points = fit.points.size();
0896 qa.z_vs_s_slope = finite_float_or_nan(fit.slope);
0897 qa.z_vs_s_intercept = finite_float_or_nan(fit.intercept);
0898 qa.z_fit_chi2 = finite_float_or_nan(fit.chi2);
0899 qa.z_fit_ndf = fit.ndf;
0900 qa.s_at_transverse_pca = finite_float_or_nan(fit.s_at_pca);
0901 qa.minimum_transverse_radius = finite_float_or_nan(fit.minimum_radius);
0902 qa.tpc_z_at_transverse_pca = finite_float_or_nan(fit.z_at_pca);
0903 qa.tpc_z_at_r0 = finite_float_or_nan(fit.z_at_r0);
0904 for (std::size_t i = 0; i < fit.points.size(); ++i)
0905 {
0906 const Point& point = fit.points[i];
0907 qa.fit_point_layer.push_back(point.layer);
0908 qa.fit_point_pad.push_back(point.pad);
0909 qa.fit_point_tbin.push_back(point.tbin);
0910 qa.fit_point_x.push_back(finite_float_or_nan(point.x));
0911 qa.fit_point_y.push_back(finite_float_or_nan(point.y));
0912 qa.fit_point_z.push_back(finite_float_or_nan(point.z));
0913 qa.fit_point_r.push_back(finite_float_or_nan(std::hypot(point.x, point.y)));
0914 qa.fit_point_s.push_back(i < fit.path_length.size() ? fit.path_length[i] : TpcCrossingCandidate::nan());
0915 }
0916
0917 std::map<unsigned int, Tpc_AssembledTrack::HitIndex> hit_by_layer;
0918 if (track)
0919 {
0920 for (unsigned int ih = 0; ih < track->size_hit_indices(); ++ih)
0921 {
0922 const Tpc_AssembledTrack::HitIndex hi = track->get_hit_index(ih);
0923 hit_by_layer.emplace(TrkrDefs::getLayer(hi.first), hi);
0924 }
0925 }
0926 if (!hit_by_layer.empty())
0927 {
0928 Point inner;
0929 if (make_xyz_point(hit_by_layer.begin()->second.first, hit_by_layer.begin()->second.second, crossing, inner))
0930 {
0931 qa.inner_layer = inner.layer;
0932 qa.inner_tbin = inner.tbin;
0933 qa.inner_x = finite_float_or_nan(inner.x);
0934 qa.inner_y = finite_float_or_nan(inner.y);
0935 qa.inner_z = finite_float_or_nan(inner.z);
0936 }
0937 Point outer;
0938 if (make_xyz_point(hit_by_layer.rbegin()->second.first, hit_by_layer.rbegin()->second.second, crossing, outer))
0939 {
0940 qa.outer_layer = outer.layer;
0941 qa.outer_tbin = outer.tbin;
0942 qa.outer_x = finite_float_or_nan(outer.x);
0943 qa.outer_y = finite_float_or_nan(outer.y);
0944 qa.outer_z = finite_float_or_nan(outer.z);
0945 }
0946 }
0947
0948 candidate.tpc_z0 = fit.z_at_r0;
0949 if (!std::isfinite(candidate.tpc_z0) || std::abs(candidate.tpc_z0 - m_collisionZ) > m_maxCandidateVertexZ)
0950 {
0951 qa.first_failed_stage = static_cast<unsigned char>(TpcCrossingCandidateStage::VertexCompatible);
0952 candidate.rejection_status = static_cast<unsigned char>(TpcCrossingStatus::VertexIncompatible);
0953 qa.rejection_status = candidate.rejection_status;
0954 return candidate;
0955 }
0956
0957 const auto vertex_iter = vertices_by_crossing.find(crossing);
0958 if (vertex_iter != vertices_by_crossing.end())
0959 {
0960 qa.n_vertices_at_crossing = vertex_iter->second.size();
0961 double best_abs_delta = std::numeric_limits<double>::max();
0962 for (const SiliconVertexHypothesis& vertex : vertex_iter->second)
0963 {
0964 const double delta_z = fit.z_at_r0 - vertex.z;
0965 const double abs_delta = std::abs(delta_z);
0966 const double pull_z = std::isfinite(vertex.sigma_z) && vertex.sigma_z > 0.0 ? delta_z / vertex.sigma_z : std::numeric_limits<double>::quiet_NaN();
0967 qa.candidate_vertex_ids.push_back(vertex.vertex_id);
0968 qa.candidate_vertex_z.push_back(finite_float_or_nan(vertex.z));
0969 qa.candidate_vertex_sigma_z.push_back(finite_float_or_nan(vertex.sigma_z));
0970 qa.candidate_vertex_ntracks.push_back(vertex.ntracks);
0971 qa.candidate_vertex_delta_z.push_back(finite_float_or_nan(delta_z));
0972 qa.candidate_vertex_pull_z.push_back(finite_float_or_nan(pull_z));
0973 if (!std::isfinite(abs_delta) || abs_delta >= best_abs_delta) { continue;
0974 }
0975
0976 best_abs_delta = abs_delta;
0977 candidate.has_silicon_vertex = true;
0978 candidate.silicon_vertex_id = vertex.vertex_id;
0979 candidate.silicon_vertex_z = vertex.z;
0980 candidate.delta_z = delta_z;
0981 qa.closest_vertex_id = vertex.vertex_id;
0982 qa.closest_vertex_z = finite_float_or_nan(vertex.z);
0983 qa.closest_vertex_sigma_z = finite_float_or_nan(vertex.sigma_z);
0984 qa.closest_vertex_delta_z = finite_float_or_nan(delta_z);
0985 qa.closest_vertex_abs_delta_z = finite_float_or_nan(abs_delta);
0986 qa.closest_vertex_pull_z = finite_float_or_nan(pull_z);
0987 }
0988 }
0989
0990 qa.has_silicon_vertex = candidate.has_silicon_vertex;
0991 if (qa.has_silicon_vertex) { qa.candidate_qa_bits |= HasSiliconVertex;
0992 }
0993 candidate.vertex_compatible = candidate.has_silicon_vertex && std::abs(candidate.delta_z) <= m_maxVertexDz;
0994 qa.vertex_compatible = candidate.vertex_compatible;
0995 if (qa.vertex_compatible) { qa.candidate_qa_bits |= PassesVertexDz;
0996 }
0997
0998 candidate.tpc_valid = true;
0999 qa.tpc_valid = candidate.tpc_valid;
1000 qa.first_failed_stage = static_cast<unsigned char>(TpcCrossingCandidateStage::FinalRanking);
1001 qa.rejection_status = candidate.rejection_status;
1002 return candidate;
1003 }
1004
1005 int TpcCrossingFinder::process_event(PHCompositeNode* topNode)
1006 {
1007 if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK) { return Fun4AllReturnCodes::ABORTEVENT;
1008 }
1009 if (!m_assembledTracks || !m_decisions) { return Fun4AllReturnCodes::EVENT_OK;
1010 }
1011 m_decisions->Reset();
1012
1013 std::set<short> available_crossings = get_available_crossings();
1014 const std::set<short> intt_crossings = get_intt_crossings();
1015 available_crossings.insert(intt_crossings.begin(), intt_crossings.end());
1016 const auto vertices_by_crossing = get_vertices_by_crossing();
1017 const bool has_vertex_map = m_vertexMap != nullptr;
1018
1019 if (Verbosity() > 0)
1020 {
1021 std::cout << Name() << "::process_event - event " << m_event << " available candidate crossings:";
1022 for (const short crossing : available_crossings) { std::cout << " " << crossing;
1023 }
1024 std::cout << " | intt crossings:";
1025 for (const short crossing : intt_crossings) { std::cout << " " << crossing;
1026 }
1027 std::cout << " | vertices by crossing:";
1028 for (const auto& item : vertices_by_crossing)
1029 {
1030 std::cout << " " << item.first << "(" << item.second.size() << ":";
1031 for (std::size_t ivtx = 0; ivtx < item.second.size(); ++ivtx)
1032 {
1033 const SiliconVertexHypothesis& vertex = item.second[ivtx];
1034 if (ivtx != 0U) { std::cout << ",";
1035 }
1036 std::cout << " id=" << vertex.vertex_id
1037 << " x=" << vertex.x
1038 << " y=" << vertex.y
1039 << " z=" << vertex.z;
1040 }
1041 std::cout << ")";
1042 }
1043 std::cout << std::endl;
1044 }
1045
1046 const unsigned int nassembled = m_assembledTracks->size();
1047 for (unsigned int iassembled = 0; iassembled < nassembled; ++iassembled)
1048 {
1049 const Tpc_AssembledTrack* assembled = m_assembledTracks->get_track(iassembled);
1050 if (!assembled) { continue;
1051 }
1052
1053 TpcCrossingDecisionv1* decision = new TpcCrossingDecisionv1();
1054 std::vector<TpcCrossingCandidate> candidate_qa_records;
1055 auto add_decision_with_candidates = [&candidate_qa_records, decision, this]()
1056 {
1057 for (const TpcCrossingCandidate& candidate : candidate_qa_records) { decision->add_candidate(candidate);
1058 }
1059 m_decisions->add_decision(decision);
1060 };
1061 decision->set_assembled_track_id(assembled->get_track_id());
1062 decision->set_number_of_available_crossings(static_cast<unsigned short>(std::min<std::size_t>(available_crossings.size(), std::numeric_limits<unsigned short>::max())));
1063
1064 if (available_crossings.empty())
1065 {
1066 decision->set_status(TpcCrossingStatus::NoInttCrossings);
1067 add_decision_with_candidates();
1068 continue;
1069 }
1070
1071 TrkrDefs::hitsetkey min_hsk = 0;
1072 TrkrDefs::hitkey min_hk = 0;
1073 TrkrDefs::hitsetkey max_hsk = 0;
1074 TrkrDefs::hitkey max_hk = 0;
1075 if (!find_time_extrema(assembled, min_hsk, min_hk, max_hsk, max_hk))
1076 {
1077 for (const short crossing : available_crossings)
1078 {
1079 TpcCrossingCandidate qa;
1080 qa.crossing = crossing;
1081 qa.is_available_from_intt = intt_crossings.contains(crossing);
1082 qa.candidate_qa_bits = qa.is_available_from_intt ? FromInttCrossing : 0U;
1083 qa.rejection_status = static_cast<unsigned char>(TpcCrossingStatus::NoValidCrossing);
1084 qa.first_failed_stage = static_cast<unsigned char>(TpcCrossingCandidateStage::PassedTimeWindow);
1085 candidate_qa_records.push_back(qa);
1086 }
1087 decision->set_status(TpcCrossingStatus::NoValidCrossing);
1088 add_decision_with_candidates();
1089 continue;
1090 }
1091
1092 const unsigned int min_tbin = TpcDefs::getTBin(min_hk);
1093 const unsigned int max_tbin = TpcDefs::getTBin(max_hk);
1094 std::vector<short> allowed_crossings;
1095 for (const short crossing : available_crossings)
1096 {
1097 const double min_time0_ns = (static_cast<double>(min_tbin) - m_t0) * m_tpcAdcClock;
1098 const double max_time0_ns = (static_cast<double>(max_tbin) - m_t0) * m_tpcAdcClock;
1099 const double min_time_ns = min_time0_ns - static_cast<double>(crossing) * m_crossingPeriodNs;
1100 const double max_time_ns = max_time0_ns - static_cast<double>(crossing) * m_crossingPeriodNs;
1101 const bool passes_time_window = std::isfinite(min_time_ns) && std::isfinite(max_time_ns) &&
1102 min_time_ns > 0.0 && max_time_ns <= m_maxLookupTimeNs;
1103 if (passes_time_window)
1104 {
1105 allowed_crossings.push_back(crossing);
1106 }
1107 else
1108 {
1109 TpcCrossingCandidate qa;
1110 qa.crossing = crossing;
1111 qa.is_available_from_intt = intt_crossings.contains(crossing);
1112 qa.rejection_status = static_cast<unsigned char>(TpcCrossingStatus::NoAllowedCrossing);
1113 qa.first_failed_stage = static_cast<unsigned char>(TpcCrossingCandidateStage::PassedTimeWindow);
1114 qa.min_tbin_time_crossing0_ns = finite_float_or_nan(min_time0_ns);
1115 qa.max_tbin_time_crossing0_ns = finite_float_or_nan(max_time0_ns);
1116 qa.candidate_min_time_ns = finite_float_or_nan(min_time_ns);
1117 qa.candidate_max_time_ns = finite_float_or_nan(max_time_ns);
1118 qa.max_lookup_time_ns = finite_float_or_nan(m_maxLookupTimeNs);
1119 qa.min_time_margin_ns = finite_float_or_nan(min_time_ns);
1120 qa.max_time_margin_ns = finite_float_or_nan(m_maxLookupTimeNs - max_time_ns);
1121 qa.min_time_layer = TrkrDefs::getLayer(min_hsk);
1122 qa.min_time_side = TpcDefs::getSide(min_hsk);
1123 qa.min_time_pad = TpcDefs::getPad(min_hk);
1124 qa.min_time_tbin = min_tbin;
1125 qa.max_time_layer = TrkrDefs::getLayer(max_hsk);
1126 qa.max_time_side = TpcDefs::getSide(max_hsk);
1127 qa.max_time_pad = TpcDefs::getPad(max_hk);
1128 qa.max_time_tbin = max_tbin;
1129 qa.candidate_qa_bits = qa.is_available_from_intt ? FromInttCrossing : 0U;
1130 candidate_qa_records.push_back(qa);
1131 }
1132 }
1133 decision->set_number_of_allowed_crossings(static_cast<unsigned short>(std::min<std::size_t>(allowed_crossings.size(), std::numeric_limits<unsigned short>::max())));
1134
1135 if (allowed_crossings.empty())
1136 {
1137 decision->set_status(TpcCrossingStatus::NoAllowedCrossing);
1138 add_decision_with_candidates();
1139 continue;
1140 }
1141
1142 std::vector<Candidate> valid_candidates;
1143 unsigned char last_rejection = static_cast<unsigned char>(TpcCrossingStatus::NoValidCrossing);
1144 for (const short crossing : allowed_crossings)
1145 {
1146 Candidate candidate = test_candidate(assembled, crossing, min_hsk, min_hk, max_hsk, max_hk, vertices_by_crossing);
1147 decision->set_number_of_tested_crossings(decision->get_number_of_tested_crossings() + 1U);
1148 candidate.qa.is_available_from_intt = intt_crossings.contains(crossing);
1149 if (!candidate.qa.is_available_from_intt) { candidate.qa.candidate_qa_bits &= ~FromInttCrossing;
1150 }
1151 candidate_qa_records.push_back(candidate.qa);
1152 if (candidate.tpc_valid)
1153 {
1154 valid_candidates.push_back(candidate);
1155 }
1156 else if (candidate.rejection_status != 0U)
1157 {
1158 last_rejection = candidate.rejection_status;
1159 }
1160 }
1161 for (Candidate& candidate : valid_candidates)
1162 {
1163 if (candidate.has_silicon_vertex)
1164 {
1165 const double abs_delta = std::isfinite(candidate.delta_z) ? std::abs(candidate.delta_z) : std::numeric_limits<double>::max();
1166 candidate.confidence_tier = abs_delta <= m_maxVertexDz ? 0U : 1U;
1167 candidate.confidence_score = abs_delta;
1168 }
1169 else
1170 {
1171 const double abs_z0 = std::isfinite(candidate.tpc_z0) ? std::abs(candidate.tpc_z0 - m_collisionZ) : std::numeric_limits<double>::max();
1172 candidate.confidence_tier = 2U;
1173 candidate.confidence_score = abs_z0;
1174 if (abs_z0 > m_maxTier2BeamlineZ)
1175 {
1176 candidate.tpc_valid = false;
1177 candidate.rejection_status = static_cast<unsigned char>(TpcCrossingStatus::VertexIncompatible);
1178 last_rejection = candidate.rejection_status;
1179 }
1180 }
1181 }
1182
1183 std::vector<Candidate> ranked_candidates;
1184 for (const Candidate& candidate : valid_candidates)
1185 {
1186 if (!candidate.tpc_valid) { continue;
1187 }
1188 if (candidate.confidence_tier == InvalidConfidenceTier) { continue;
1189 }
1190 if (!std::isfinite(candidate.confidence_score)) { continue;
1191 }
1192 ranked_candidates.push_back(candidate);
1193 }
1194
1195 for (TpcCrossingCandidate& qa : candidate_qa_records)
1196 {
1197 const auto iq = std::find_if(valid_candidates.begin(), valid_candidates.end(), [&qa](const Candidate& candidate) { return candidate.crossing == qa.crossing; });
1198 if (iq == valid_candidates.end()) { continue;
1199 }
1200 qa.confidence_tier = iq->confidence_tier;
1201 qa.confidence_score = finite_float_or_nan(iq->confidence_score);
1202 qa.tpc_valid = iq->tpc_valid;
1203 if (!iq->tpc_valid && iq->rejection_status != 0U) { qa.rejection_status = iq->rejection_status;
1204 }
1205 }
1206
1207 decision->set_number_of_tpc_valid_crossings(static_cast<unsigned short>(std::min<std::size_t>(ranked_candidates.size(), std::numeric_limits<unsigned short>::max())));
1208 const auto n_vertex_compatible = std::count_if(ranked_candidates.begin(), ranked_candidates.end(), [](const Candidate& candidate) { return candidate.confidence_tier == 0U; });
1209 decision->set_number_of_vertex_compatible_crossings(static_cast<unsigned short>(std::min<std::size_t>(n_vertex_compatible, std::numeric_limits<unsigned short>::max())));
1210
1211 if (ranked_candidates.empty())
1212 {
1213 decision->set_status(last_rejection != 0U ? last_rejection : static_cast<unsigned char>(TpcCrossingStatus::NoValidCrossing));
1214 add_decision_with_candidates();
1215 continue;
1216 }
1217
1218 std::sort(ranked_candidates.begin(), ranked_candidates.end(), [](const Candidate& a, const Candidate& b) {
1219 if (a.confidence_tier != b.confidence_tier) { return a.confidence_tier < b.confidence_tier;
1220 }
1221 if (a.confidence_score != b.confidence_score) { return a.confidence_score < b.confidence_score;
1222 }
1223 return a.crossing < b.crossing;
1224 });
1225
1226 const Candidate& best_candidate = ranked_candidates.front();
1227 const double best_score = best_candidate.confidence_score;
1228 double second_score = std::numeric_limits<double>::infinity();
1229 for (std::size_t icandidate = 1; icandidate < ranked_candidates.size(); ++icandidate)
1230 {
1231 if (ranked_candidates[icandidate].confidence_tier != best_candidate.confidence_tier) { break;
1232 }
1233 second_score = ranked_candidates[icandidate].confidence_score;
1234 break;
1235 }
1236 const bool close_rival = std::isfinite(second_score) && second_score - best_score <= m_minBestSecondSeparation;
1237
1238 auto selected_status = TpcCrossingStatus::SelectedByContainment;
1239 if (best_candidate.confidence_tier == 0U)
1240 {
1241 selected_status = close_rival ? TpcCrossingStatus::SelectedByVertexAmbiguous : TpcCrossingStatus::SelectedByVertex;
1242 }
1243 else if (best_candidate.confidence_tier == 1U)
1244 {
1245 selected_status = TpcCrossingStatus::SelectedByVertexLoose;
1246 }
1247 else if (best_candidate.confidence_tier == 2U)
1248 {
1249 selected_status = close_rival ? TpcCrossingStatus::SelectedByContainmentAmbiguous : TpcCrossingStatus::SelectedByContainment;
1250 }
1251
1252 decision->set_best_abs_delta_z(finite_float_or_nan(best_score));
1253 decision->set_second_best_abs_delta_z(finite_float_or_nan(second_score));
1254 decision->set_selected_crossing(best_candidate.crossing);
1255 decision->set_tpc_z0(finite_float_or_nan(best_candidate.tpc_z0));
1256 decision->set_silicon_vertex_id(best_candidate.silicon_vertex_id);
1257 decision->set_silicon_vertex_z(finite_float_or_nan(best_candidate.silicon_vertex_z));
1258 decision->set_delta_z(finite_float_or_nan(best_candidate.delta_z));
1259 decision->set_selected_tier(best_candidate.confidence_tier);
1260 decision->set_selected_score(finite_float_or_nan(best_candidate.confidence_score));
1261 decision->set_status(selected_status);
1262
1263 std::vector<short> by_vertex;
1264 by_vertex.reserve(valid_candidates.size());
1265 std::vector<short> by_collision;
1266 by_collision.reserve(valid_candidates.size());
1267 for (const Candidate& candidate : valid_candidates)
1268 {
1269 by_vertex.push_back(candidate.crossing);
1270 by_collision.push_back(candidate.crossing);
1271 }
1272 std::sort(by_vertex.begin(), by_vertex.end(), [&valid_candidates](short a, short b) {
1273 const auto ia = std::find_if(valid_candidates.begin(), valid_candidates.end(), [a](const Candidate& c) { return c.crossing == a; });
1274 const auto ib = std::find_if(valid_candidates.begin(), valid_candidates.end(), [b](const Candidate& c) { return c.crossing == b; });
1275 const double da = std::isfinite(ia->delta_z) ? std::abs(ia->delta_z) : std::numeric_limits<double>::max();
1276 const double db = std::isfinite(ib->delta_z) ? std::abs(ib->delta_z) : std::numeric_limits<double>::max();
1277 if (da != db) { return da < db;
1278 }
1279 return a < b;
1280 });
1281 std::sort(by_collision.begin(), by_collision.end(), [this, &valid_candidates](short a, short b) {
1282 const auto ia = std::find_if(valid_candidates.begin(), valid_candidates.end(), [a](const Candidate& c) { return c.crossing == a; });
1283 const auto ib = std::find_if(valid_candidates.begin(), valid_candidates.end(), [b](const Candidate& c) { return c.crossing == b; });
1284 const double da = std::isfinite(ia->tpc_z0) ? std::abs(ia->tpc_z0 - m_collisionZ) : std::numeric_limits<double>::max();
1285 const double db = std::isfinite(ib->tpc_z0) ? std::abs(ib->tpc_z0 - m_collisionZ) : std::numeric_limits<double>::max();
1286 if (da != db) { return da < db;
1287 }
1288 return a < b;
1289 });
1290
1291 for (TpcCrossingCandidate& qa : candidate_qa_records)
1292 {
1293 const auto iv = std::find(by_vertex.begin(), by_vertex.end(), qa.crossing);
1294 if (iv != by_vertex.end()) { qa.candidate_rank_by_abs_delta_z = std::distance(by_vertex.begin(), iv);
1295 }
1296 const auto ic = std::find(by_collision.begin(), by_collision.end(), qa.crossing);
1297 if (ic != by_collision.end()) { qa.candidate_rank_by_abs_collision_z = std::distance(by_collision.begin(), ic);
1298 }
1299 const auto iq = std::find_if(valid_candidates.begin(), valid_candidates.end(), [&qa](const Candidate& candidate) { return candidate.crossing == qa.crossing; });
1300 if (iq != valid_candidates.end())
1301 {
1302 qa.confidence_tier = iq->confidence_tier;
1303 qa.confidence_score = finite_float_or_nan(iq->confidence_score);
1304 qa.tpc_valid = iq->tpc_valid;
1305 if (!iq->tpc_valid && iq->rejection_status != 0U) { qa.rejection_status = iq->rejection_status;
1306 }
1307 }
1308 qa.is_selected = qa.crossing == decision->get_selected_crossing();
1309 if (qa.is_selected)
1310 {
1311 qa.candidate_qa_bits |= IsSelected;
1312 qa.first_failed_stage = static_cast<unsigned char>(TpcCrossingCandidateStage::Selected);
1313 }
1314 if (qa.candidate_rank_by_abs_delta_z == 0) { qa.candidate_qa_bits |= IsBestVertexCandidate;
1315 }
1316 }
1317
1318 add_decision_with_candidates();
1319 }
1320
1321 if (Verbosity() > 0)
1322 {
1323 std::cout << Name() << "::process_event - event " << m_event
1324 << " assembled_tracks=" << nassembled
1325 << " decisions=" << m_decisions->size()
1326 << " available_crossings=" << available_crossings.size()
1327 << " vertex_map=" << has_vertex_map
1328 << std::endl;
1329 }
1330
1331 ++m_event;
1332 return Fun4AllReturnCodes::EVENT_OK;
1333 }
