File indexing completed on 2026-08-31 08:21:23
0001 #include "Tpc_PolyClusterizer.h"
0002
0003 #include "IdealPadMap.h"
0004 #include "Tpc_AssembledTrack.h"
0005 #include "Tpc_AssembledTrackContainer.h"
0006 #include "Tpc_PolyClusterContainerv1.h"
0007 #include "Tpc_PolyClusterv1.h"
0008 #include "TpcCrossingDecision.h"
0009 #include "TpcCrossingDecisionContainer.h"
0010
0011 #include <cdbobjects/CDBTTree.h>
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013
0014 #include <ffamodules/CDBInterface.h>
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 <trackbase/ActsGeometry.h>
0028 #include <g4detectors/PHG4CylinderGeom.h> // for PHG4CylinderGeom
0029 #include <g4detectors/PHG4CylinderGeomContainer.h>
0030 #include <g4detectors/PHG4TpcGeom.h>
0031 #include <g4detectors/PHG4TpcGeomContainer.h>
0032
0033
0034 #include <TPolyLine3D.h>
0035 #include <phgarfield/PHGarfield.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 <vector>
0045
0046 namespace
0047 {
0048 constexpr unsigned int FirstLayer = 7;
0049 constexpr unsigned int LastLayer = 54;
0050 constexpr unsigned int NLayers = LastLayer - FirstLayer + 1;
0051 constexpr unsigned int NSides = 2;
0052 constexpr unsigned int NSectors = 12;
0053 constexpr double PhiConsistencyTolerance = 1.0e-10;
0054
0055 double wrap_phi(double phi)
0056 {
0057 while (phi > M_PI)
0058 {
0059 phi -= 2.0 * M_PI;
0060 }
0061 while (phi <= -M_PI)
0062 {
0063 phi += 2.0 * M_PI;
0064 }
0065 return phi;
0066 }
0067
0068 double unwrap_phi_near(double phi, const double reference)
0069 {
0070 while (phi - reference > M_PI)
0071 {
0072 phi -= 2.0 * M_PI;
0073 }
0074 while (phi - reference <= -M_PI)
0075 {
0076 phi += 2.0 * M_PI;
0077 }
0078 return phi;
0079 }
0080
0081 double clamp_unit(const double value)
0082 {
0083 return std::max(0.0, std::min(1.0, value));
0084 }
0085
0086 double square(const double value)
0087 {
0088 return value * value;
0089 }
0090
0091 double phi_sample_fraction(const unsigned int sample)
0092 {
0093 if (Tpc_PolyClusterizer::NPhiSamples <= 1U)
0094 {
0095 return 0.0;
0096 }
0097 return static_cast<double>(sample) / static_cast<double>(Tpc_PolyClusterizer::NPhiSamples - 1U);
0098 }
0099 }
0100
0101 Tpc_PolyClusterizer::Tpc_PolyClusterizer(const std::string& name)
0102 : SubsysReco(name)
0103 , m_inputNodeName("TPC_ASSEMBLEDTRACKS")
0104 , m_outputNodeName("TPC_POLYCLUSTERS")
0105 {
0106 }
0107
0108 Tpc_PolyClusterizer::~Tpc_PolyClusterizer()
0109 {
0110 delete m_idealPadMap;
0111 m_idealPadMap = nullptr;
0112 delete m_garfield;
0113 m_garfield = nullptr;
0114 }
0115
0116 int Tpc_PolyClusterizer::InitRun(PHCompositeNode* topNode)
0117 {
0118 if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0119 {
0120 return Fun4AllReturnCodes::ABORTRUN;
0121 }
0122 if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0123 {
0124 return Fun4AllReturnCodes::ABORTRUN;
0125 }
0126
0127 delete m_idealPadMap;
0128 m_idealPadMap = new IdealPadMap();
0129 if (m_idealPadMap->load_from_cdb(Verbosity()) != 0 || !m_idealPadMap->is_loaded())
0130 {
0131 std::cerr << Name() << "::InitRun - failed to load IdealPadMap" << std::endl;
0132 return Fun4AllReturnCodes::ABORTRUN;
0133 }
0134
0135 PHG4TpcGeom *layergeom = m_geomContainerTpc->GetLayerCellGeom(20);
0136 double rot_x = layergeom->get_rot_x();
0137 double rot_y = layergeom->get_rot_y();
0138 double rot_z = layergeom->get_rot_z();
0139 double place_x = layergeom->get_place_x();
0140 double place_y = layergeom->get_place_y();
0141 double place_z = layergeom->get_place_z();
0142 if (use_survey_geometry)
0143 {
0144 m_tpcMove = {place_x, place_y, place_z};
0145 m_tpcRotations = {{{rot_x, rot_y, rot_z}, {0.0, 0.0, 0.0}}};
0146 }
0147
0148 delete m_garfield;
0149
0150 const std::string electricFieldMap = CDBInterface::instance()->getUrl("Tpc_PolySeeding_EField");
0151
0152
0153
0154 const auto kefffile = CDBInterface::instance()->getUrl("Tpc_PolyClusterizer_kEff");
0155
0156 if (!kefffile.empty())
0157 {
0158 auto keffcdbtree = std::make_unique<CDBTTree>(kefffile);
0159 keffcdbtree->LoadCalibrations();
0160 m_kEffSide0 = keffcdbtree->GetSingleFloatValue("keffside0");
0161 m_kEffSide1 = keffcdbtree->GetSingleFloatValue("keffside1");
0162
0163 std::cout << PHWHERE << "Poly clusterizer loading calibrations from " << kefffile << std::endl
0164 <<" with values keffside0 = " << m_kEffSide0 << ", keffside1 = " << m_kEffSide1 << std::endl;
0165 }
0166
0167
0168 m_garfield = new PHGarfield(Name() + "_PHGarfield", electricFieldMap, m_kEffSide0, m_kEffSide1);
0169 configure_garfield(m_garfield);
0170 if (m_garfield->InitRun(topNode) != Fun4AllReturnCodes::EVENT_OK)
0171 {
0172 std::cerr << Name() << "::InitRun - PHGarfield InitRun failed" << std::endl;
0173 return Fun4AllReturnCodes::ABORTRUN;
0174 }
0175
0176 if (!build_drift_lookup())
0177 {
0178 return Fun4AllReturnCodes::ABORTRUN;
0179 }
0180
0181 m_event = 0;
0182 return Fun4AllReturnCodes::EVENT_OK;
0183 }
0184
0185 void Tpc_PolyClusterizer::configure_garfield(PHGarfield* garfield) const
0186 {
0187 if (!garfield)
0188 {
0189 return;
0190 }
0191
0192 garfield->MoveTpc(m_tpcMove[0], m_tpcMove[1], m_tpcMove[2]);
0193 for (const auto& rotation : m_tpcRotations)
0194 {
0195 garfield->RotateTpc(rotation[0], rotation[1], rotation[2]);
0196 }
0197 garfield->SetCMVoltageDefault(m_cmVoltageDefault);
0198 }
0199
0200 int Tpc_PolyClusterizer::getNodes(PHCompositeNode* topNode)
0201 {
0202 m_assembledTracks = findNode::getClass<Tpc_AssembledTrackContainer>(topNode, m_inputNodeName);
0203 if (!m_assembledTracks)
0204 {
0205 std::cerr << Name() << "::getNodes - missing " << m_inputNodeName << std::endl;
0206 return Fun4AllReturnCodes::ABORTRUN;
0207 }
0208
0209 m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0210 if (!m_hits)
0211 {
0212 std::cerr << Name() << "::getNodes - missing TRKR_HITSET" << std::endl;
0213 return Fun4AllReturnCodes::ABORTRUN;
0214 }
0215
0216 m_crossingDecisions = findNode::getClass<TpcCrossingDecisionContainer>(topNode, m_crossingDecisionNodeName);
0217 if (!m_crossingDecisions)
0218 {
0219 std::cerr << Name() << "::getNodes - missing " << m_crossingDecisionNodeName << std::endl;
0220 return Fun4AllReturnCodes::ABORTRUN;
0221 }
0222
0223 m_geomContainerTpc = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0224 if (!m_geomContainerTpc)
0225 {
0226 std::cerr << Name() << "::getNodes - missing TPCGEOMCONTAINER" << std::endl;
0227 return Fun4AllReturnCodes::ABORTRUN;
0228 }
0229
0230 return Fun4AllReturnCodes::EVENT_OK;
0231 }
0232
0233 int Tpc_PolyClusterizer::createNodes(PHCompositeNode* topNode)
0234 {
0235 PHNodeIterator iter(topNode);
0236 PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0237 if (!dstNode)
0238 {
0239 dstNode = new PHCompositeNode("DST");
0240 topNode->addNode(dstNode);
0241 }
0242
0243 m_clusters = findNode::getClass<Tpc_PolyClusterContainer>(topNode, m_outputNodeName);
0244 if (!m_clusters)
0245 {
0246 m_clusters = new Tpc_PolyClusterContainerv1();
0247 PHIODataNode<PHObject>* node = new PHIODataNode<PHObject>(m_clusters, m_outputNodeName, "PHObject");
0248 dstNode->addNode(node);
0249 std::cout << Name() << "::createNodes - created " << m_outputNodeName << " node" << std::endl;
0250 }
0251
0252 return Fun4AllReturnCodes::EVENT_OK;
0253 }
0254
0255 unsigned int Tpc_PolyClusterizer::drift_lookup_index(const unsigned int layer_index,
0256 const unsigned int side,
0257 const unsigned int sector,
0258 const unsigned int sample)
0259 {
0260 return (((layer_index * NSides + side) * NSectors + sector) * NPhiSamples + sample);
0261 }
0262
0263 bool Tpc_PolyClusterizer::build_drift_lookup()
0264 {
0265 if (!m_idealPadMap || !m_garfield)
0266 {
0267 return false;
0268 }
0269 if (m_reverseDriftStepNs <= 0.0 || !std::isfinite(m_reverseDriftStepNs))
0270 {
0271 return false;
0272 }
0273
0274 for (DriftPolyline& polyline : m_driftLookup)
0275 {
0276 polyline.phi = 0.0;
0277 polyline.points.clear();
0278 }
0279
0280 unsigned int nbuilt = 0;
0281 for (unsigned int layer = FirstLayer; layer <= LastLayer; ++layer)
0282 {
0283 const unsigned int layer_index = layer - FirstLayer;
0284 const double radius = m_idealPadMap->get_radius(layer);
0285 const unsigned int pads_per_sector = m_idealPadMap->get_pads_per_sector_for_layer(layer);
0286 if (!std::isfinite(radius) || pads_per_sector == 0U)
0287 {
0288 std::cerr << Name() << "::build_drift_lookup - invalid geometry for layer " << layer << std::endl;
0289 return false;
0290 }
0291
0292 for (unsigned int side = 0; side < NSides; ++side)
0293 {
0294 const double z0 = (side == 0U) ? m_startZSouth : m_startZNorth;
0295 for (unsigned int sector = 0; sector < NSectors; ++sector)
0296 {
0297 for (unsigned int sample = 0; sample < NPhiSamples; ++sample)
0298 {
0299 const unsigned int local_phibin = static_cast<unsigned int>(std::llround(
0300 phi_sample_fraction(sample) * static_cast<double>(pads_per_sector - 1U)));
0301 const unsigned int global_pad = sector * pads_per_sector + local_phibin;
0302 const double phi_local = m_idealPadMap->get_phi(side, sector, layer, local_phibin);
0303 const double phi_global = m_idealPadMap->get_phi(side, layer, global_pad);
0304 if (!std::isfinite(phi_local) || !std::isfinite(phi_global))
0305 {
0306 std::cerr << Name() << "::build_drift_lookup - invalid phi for layer " << layer
0307 << " side " << side
0308 << " sector " << sector
0309 << " local_phibin " << local_phibin
0310 << " global_pad " << global_pad
0311 << " phi_local " << phi_local
0312 << " phi_global " << phi_global << std::endl;
0313 return false;
0314 }
0315
0316 const double phi_difference = wrap_phi(phi_global - phi_local);
0317 if (std::abs(phi_difference) > PhiConsistencyTolerance)
0318 {
0319 std::cerr << Name() << "::build_drift_lookup - inconsistent IdealPadMap phi overloads"
0320 << " layer " << layer
0321 << " side " << side
0322 << " sector " << sector
0323 << " local_phibin " << local_phibin
0324 << " global_pad " << global_pad
0325 << " phi_local " << phi_local
0326 << " phi_global " << phi_global
0327 << " phi_difference " << phi_difference << std::endl;
0328 return false;
0329 }
0330
0331 const double phi = phi_local;
0332 const double x0 = radius * std::cos(phi);
0333 const double y0 = radius * std::sin(phi);
0334 TPolyLine3D* drift = m_garfield->ReverseDrift(x0, y0, z0, m_reverseDriftStepNs);
0335 if (!drift || drift->GetN() <= 0)
0336 {
0337 delete drift;
0338 std::cerr << Name() << "::build_drift_lookup - ReverseDrift failed for layer " << layer
0339 << " side " << side << " sector " << sector << " sample " << sample << std::endl;
0340 return false;
0341 }
0342
0343 const int npoints = drift->GetN();
0344 const Float_t* xyz = drift->GetP();
0345 if (!xyz || npoints <= 0)
0346 {
0347 delete drift;
0348 std::cerr << Name() << "::build_drift_lookup - empty drift points for layer " << layer
0349 << " side " << side << " sector " << sector << " sample " << sample << std::endl;
0350 return false;
0351 }
0352
0353 DriftPolyline& polyline = m_driftLookup[drift_lookup_index(layer_index, side, sector, sample)];
0354 polyline.phi = phi;
0355 polyline.points.resize(static_cast<std::size_t>(npoints));
0356 for (int ipoint = 0; ipoint < npoints; ++ipoint)
0357 {
0358 const int idx = 3 * ipoint;
0359 DriftPoint& point = polyline.points[static_cast<std::size_t>(ipoint)];
0360 const double output_r = std::hypot(static_cast<double>(xyz[idx]), static_cast<double>(xyz[idx + 1]));
0361 const double output_phi = unwrap_phi_near(std::atan2(static_cast<double>(xyz[idx + 1]), static_cast<double>(xyz[idx])), phi);
0362 point.delta_r = static_cast<float>(output_r - radius);
0363 point.delta_phi = static_cast<float>(output_phi - phi);
0364 point.z = xyz[idx + 2];
0365 }
0366 delete drift;
0367 ++nbuilt;
0368 }
0369 }
0370 }
0371 }
0372
0373 for (unsigned int layer = FirstLayer; layer <= LastLayer; ++layer)
0374 {
0375 const unsigned int layer_index = layer - FirstLayer;
0376 const double radius = m_idealPadMap->get_radius(layer);
0377 for (unsigned int side = 0; side < NSides; ++side)
0378 {
0379 const DriftPolyline& reference = m_driftLookup[drift_lookup_index(layer_index, side, 0, 0)];
0380 bool symmetry_ok = !reference.points.empty() && std::isfinite(radius);
0381 unsigned int length_mismatches = 0;
0382 double max_delta_r = 0.0;
0383 double max_delta_phi_arc = 0.0;
0384 double max_delta_z = 0.0;
0385
0386 for (unsigned int sector = 0; sector < NSectors; ++sector)
0387 {
0388 for (unsigned int sample = 0; sample < NPhiSamples; ++sample)
0389 {
0390 const DriftPolyline& polyline = m_driftLookup[drift_lookup_index(layer_index, side, sector, sample)];
0391 if (polyline.points.size() != reference.points.size())
0392 {
0393 ++length_mismatches;
0394 symmetry_ok = false;
0395 }
0396
0397 const std::size_t npoints = std::min(reference.points.size(), polyline.points.size());
0398 for (std::size_t ipoint = 0; ipoint < npoints; ++ipoint)
0399 {
0400 const DriftPoint& ref_point = reference.points[ipoint];
0401 const DriftPoint& point = polyline.points[ipoint];
0402 max_delta_r = std::max(max_delta_r, std::abs(static_cast<double>(point.delta_r - ref_point.delta_r)));
0403 max_delta_phi_arc = std::max(max_delta_phi_arc,
0404 std::abs(radius * wrap_phi(static_cast<double>(point.delta_phi - ref_point.delta_phi))));
0405 max_delta_z = std::max(max_delta_z, std::abs(static_cast<double>(point.z - ref_point.z)));
0406 }
0407 }
0408 }
0409
0410 if (max_delta_r > 0.001 ||
0411 max_delta_phi_arc > 0.001 ||
0412 max_delta_z > 0.001)
0413 {
0414 symmetry_ok = false;
0415 }
0416
0417 std::cout << Name() << "::build_drift_lookup - phi symmetry "
0418 << (symmetry_ok ? "holds" : "broken")
0419 << " layer=" << layer
0420 << " side=" << side
0421 << " max_dr=" << max_delta_r
0422 << " max_r_dphi=" << max_delta_phi_arc
0423 << " max_dz=" << max_delta_z
0424 << " length_mismatches=" << length_mismatches << std::endl;
0425 }
0426 }
0427
0428 if (Verbosity() > 0)
0429 {
0430 std::cout << Name() << "::build_drift_lookup - built " << nbuilt << " drift polylines" << std::endl;
0431 }
0432 return nbuilt == NLayers * NSides * NSectors * NPhiSamples;
0433 }
0434
0435 bool Tpc_PolyClusterizer::sample_drift_lookup(const unsigned int layer,
0436 const unsigned int side,
0437 const unsigned int pad,
0438 const unsigned int tbin,
0439 const short crossing,
0440 double& x,
0441 double& y,
0442 double& z) const
0443 {
0444 if (!m_idealPadMap)
0445 {
0446 return false;
0447 }
0448 if (layer < FirstLayer || layer > LastLayer)
0449 {
0450 return false;
0451 }
0452 if (side >= NSides)
0453 {
0454 return false;
0455 }
0456 if (m_reverseDriftStepNs <= 0.0 || !std::isfinite(m_reverseDriftStepNs))
0457 {
0458 return false;
0459 }
0460
0461 const unsigned int pads_per_sector = m_idealPadMap->get_pads_per_sector_for_layer(layer);
0462 if (pads_per_sector == 0U)
0463 {
0464 return false;
0465 }
0466
0467 const unsigned int sector = pad / pads_per_sector;
0468 if (sector >= NSectors)
0469 {
0470 return false;
0471 }
0472
0473 const double hit_radius = m_idealPadMap->get_radius(layer);
0474 const double hit_phi = m_idealPadMap->get_phi(side, layer, pad);
0475 if (!std::isfinite(hit_radius) || !std::isfinite(hit_phi))
0476 {
0477 return false;
0478 }
0479
0480 const double target_time_ns = (static_cast<double>(tbin) - m_t0) * m_tpcAdcClock - static_cast<double>(crossing) * m_crossingPeriodNs;
0481 if (target_time_ns <= 0.0 || !std::isfinite(target_time_ns))
0482 {
0483 return false;
0484 }
0485
0486 const unsigned int layer_index = layer - FirstLayer;
0487 std::array<const DriftPolyline*, NPhiSamples> samples{};
0488 std::array<double, NPhiSamples> sample_phi{};
0489 for (unsigned int sample = 0; sample < NPhiSamples; ++sample)
0490 {
0491 samples[sample] = &m_driftLookup[drift_lookup_index(layer_index, side, sector, sample)];
0492 if (samples[sample]->points.empty())
0493 {
0494 return false;
0495 }
0496 sample_phi[sample] = samples[sample]->phi;
0497 if (sample > 0U)
0498 {
0499 sample_phi[sample] = unwrap_phi_near(sample_phi[sample], sample_phi[sample - 1U]);
0500 }
0501 }
0502
0503 const double unwrapped_hit_phi = unwrap_phi_near(hit_phi, sample_phi[NPhiSamples / 2U]);
0504 const bool increasing = sample_phi[NPhiSamples - 1U] >= sample_phi[0];
0505
0506 bool bracket_found = false;
0507 unsigned int sample0 = 0;
0508 unsigned int sample1 = 0;
0509 double phi_fraction = 0.0;
0510 if ((increasing && unwrapped_hit_phi <= sample_phi[0]) || (!increasing && unwrapped_hit_phi >= sample_phi[0]))
0511 {
0512 sample0 = 0;
0513 sample1 = 0;
0514 bracket_found = true;
0515 }
0516 else if ((increasing && unwrapped_hit_phi >= sample_phi[NPhiSamples - 1U]) ||
0517 (!increasing && unwrapped_hit_phi <= sample_phi[NPhiSamples - 1U]))
0518 {
0519 sample0 = NPhiSamples - 1U;
0520 sample1 = NPhiSamples - 1U;
0521 bracket_found = true;
0522 }
0523 else
0524 {
0525 for (unsigned int sample = 0; sample + 1U < NPhiSamples; ++sample)
0526 {
0527 const bool in_interval = increasing ? (unwrapped_hit_phi >= sample_phi[sample] && unwrapped_hit_phi <= sample_phi[sample + 1U]) : (unwrapped_hit_phi <= sample_phi[sample] && unwrapped_hit_phi >= sample_phi[sample + 1U]);
0528 if (!in_interval)
0529 {
0530 continue;
0531 }
0532
0533 sample0 = sample;
0534 sample1 = sample + 1U;
0535 const double denom = sample_phi[sample1] - sample_phi[sample0];
0536 phi_fraction = (denom != 0.0) ? clamp_unit((unwrapped_hit_phi - sample_phi[sample0]) / denom) : 0.0;
0537 bracket_found = true;
0538 break;
0539 }
0540 }
0541
0542 if (!bracket_found)
0543 {
0544 std::cerr << Name() << "::sample_drift_lookup - failed to bracket hit phi"
0545 << " layer " << layer
0546 << " side " << side
0547 << " sector " << sector
0548 << " pad " << pad
0549 << " hit_phi " << hit_phi
0550 << " unwrapped_hit_phi " << unwrapped_hit_phi
0551 << " sample_phi";
0552 for (unsigned int sample = 0; sample < NPhiSamples; ++sample)
0553 {
0554 std::cerr << " " << sample_phi[sample];
0555 }
0556 std::cerr << std::endl;
0557 return false;
0558 }
0559
0560 auto sample_time = [this, target_time_ns](const DriftPolyline& polyline,
0561 double& delta_r,
0562 double& delta_phi,
0563 double& point_z) -> bool
0564 {
0565 const int npoints = static_cast<int>(polyline.points.size());
0566 if (npoints <= 0)
0567 {
0568 return false;
0569 }
0570
0571 const double max_time_ns = static_cast<double>(npoints - 1) * m_reverseDriftStepNs;
0572 if (target_time_ns > max_time_ns)
0573 {
0574 return false;
0575 }
0576
0577 const double fbin = target_time_ns / m_reverseDriftStepNs;
0578 const int i0 = std::min(static_cast<int>(std::floor(fbin)), npoints - 1);
0579 const int i1 = std::min(i0 + 1, npoints - 1);
0580 const double frac = fbin - static_cast<double>(i0);
0581
0582 const DriftPoint& p0 = polyline.points[static_cast<std::size_t>(i0)];
0583 const DriftPoint& p1 = polyline.points[static_cast<std::size_t>(i1)];
0584 const double dphi0 = static_cast<double>(p0.delta_phi);
0585 const double dphi1 = unwrap_phi_near(static_cast<double>(p1.delta_phi), dphi0);
0586 delta_r = static_cast<double>(p0.delta_r) + frac * static_cast<double>(p1.delta_r - p0.delta_r);
0587 delta_phi = dphi0 + frac * (dphi1 - dphi0);
0588 point_z = static_cast<double>(p0.z) + frac * static_cast<double>(p1.z - p0.z);
0589 return std::isfinite(delta_r) && std::isfinite(delta_phi) && std::isfinite(point_z);
0590 };
0591
0592 double delta_r0 = 0.0;
0593 double delta_phi0 = 0.0;
0594 double z0 = 0.0;
0595 const bool valid0 = sample_time(*samples[sample0], delta_r0, delta_phi0, z0);
0596
0597 double delta_r1 = 0.0;
0598 double delta_phi1 = 0.0;
0599 double z1 = 0.0;
0600 const bool same_sample = sample0 == sample1;
0601 const bool valid1 = same_sample ? valid0 : sample_time(*samples[sample1], delta_r1, delta_phi1, z1);
0602 if (same_sample)
0603 {
0604 delta_r1 = delta_r0;
0605 delta_phi1 = delta_phi0;
0606 z1 = z0;
0607 }
0608
0609 if (!valid0 && !valid1)
0610 {
0611 return false;
0612 }
0613
0614 double delta_r = 0.0;
0615 double delta_phi = 0.0;
0616 if (valid0 && valid1 && !same_sample)
0617 {
0618 delta_r = delta_r0 + phi_fraction * (delta_r1 - delta_r0);
0619 const double unwrapped_delta_phi1 = unwrap_phi_near(delta_phi1, delta_phi0);
0620 delta_phi = delta_phi0 + phi_fraction * (unwrapped_delta_phi1 - delta_phi0);
0621 z = z0 + phi_fraction * (z1 - z0);
0622 }
0623 else if (valid0)
0624 {
0625 delta_r = delta_r0;
0626 delta_phi = delta_phi0;
0627 z = z0;
0628 }
0629 else
0630 {
0631 delta_r = delta_r1;
0632 delta_phi = delta_phi1;
0633 z = z1;
0634 }
0635
0636 const double radius = hit_radius + delta_r;
0637 const double output_phi = unwrapped_hit_phi + delta_phi;
0638 x = radius * std::cos(output_phi);
0639 y = radius * std::sin(output_phi);
0640 return std::isfinite(x) && std::isfinite(y) && std::isfinite(z);
0641 }
0642
0643 bool Tpc_PolyClusterizer::make_xyz_point(TrkrDefs::hitsetkey hsk,
0644 TrkrDefs::hitkey hk,
0645 const short crossing,
0646 Point& p) const
0647 {
0648 if (!m_hits || !m_idealPadMap)
0649 {
0650 return false;
0651 }
0652
0653 TrkrHitSet* hitset = m_hits->findHitSet(hsk);
0654 if (!hitset)
0655 {
0656 return false;
0657 }
0658 TrkrHit* hit = hitset->getHit(hk);
0659 if (!hit)
0660 {
0661 return false;
0662 }
0663
0664 const unsigned int layer = TrkrDefs::getLayer(hsk);
0665 const unsigned int hit_side = TpcDefs::getSide(hsk);
0666 const unsigned int pad = TpcDefs::getPad(hk);
0667 const unsigned int tbin = TpcDefs::getTBin(hk);
0668 const double adc = hit->getAdc();
0669
0670 double x = 0.0;
0671 double y = 0.0;
0672 double z = 0.0;
0673 if (!sample_drift_lookup(layer, hit_side, pad, tbin, crossing, x, y, z))
0674 {
0675 return false;
0676 }
0677
0678 p.hitsetkey = hsk;
0679 p.hitkey = hk;
0680 p.layer = layer;
0681 p.side = hit_side;
0682 p.pad = pad;
0683 p.tbin = tbin;
0684 p.adc = adc;
0685 p.x = x;
0686 p.y = y;
0687 p.z = z;
0688 return true;
0689 }
0690
0691 Tpc_PolyClusterizer::ClusterParameters
0692 Tpc_PolyClusterizer::make_cluster_parameters(const std::vector<Point>& points,
0693 const Centroid& centroid,
0694 const int side) const
0695 {
0696 ClusterParameters params;
0697 if (points.empty() || !centroid.ok || !m_idealPadMap)
0698 {
0699 return params;
0700 }
0701
0702 std::set<unsigned int> pads;
0703 std::set<unsigned int> tbins;
0704 std::map<unsigned int, double> adc_by_pad;
0705 for (const Point& p : points)
0706 {
0707 params.adc += p.adc;
0708 pads.insert(p.pad);
0709 tbins.insert(p.tbin);
0710 adc_by_pad[p.pad] += p.adc;
0711 }
0712
0713 params.phi_width = static_cast<unsigned int>(pads.size());
0714 params.time_width = static_cast<unsigned int>(tbins.size());
0715
0716 unsigned int max_adc_pad = 0;
0717 double max_adc = -std::numeric_limits<double>::max();
0718 for (const auto& pad_adc : adc_by_pad)
0719 {
0720 if (pad_adc.second > max_adc)
0721 {
0722 max_adc = pad_adc.second;
0723 max_adc_pad = pad_adc.first;
0724 }
0725 }
0726
0727 const unsigned int total_phibins = m_idealPadMap->get_total_phibins(centroid.layer);
0728 const double pad_phi_width = total_phibins > 0U ? 2.0 * M_PI / static_cast<double>(total_phibins) : 0.0;
0729 const double cluster_phi = std::atan2(centroid.y, centroid.x);
0730 const double max_adc_phi = m_idealPadMap->get_phi(static_cast<unsigned int>(side), centroid.layer, max_adc_pad);
0731 if (pad_phi_width > 0.0 && std::isfinite(cluster_phi) && std::isfinite(max_adc_phi))
0732 {
0733 params.phase = wrap_phi(cluster_phi - max_adc_phi) / pad_phi_width;
0734 }
0735
0736 return params;
0737 }
0738
0739 Tpc_PolyClusterizer::Centroid
0740 Tpc_PolyClusterizer::make_centroid(const std::vector<Point>& points)
0741 {
0742 Centroid c;
0743 if (points.empty())
0744 {
0745 return c;
0746 }
0747
0748 double sx = 0.0;
0749 double sy = 0.0;
0750 double sz = 0.0;
0751 for (const Point& p : points)
0752 {
0753 sx += p.x;
0754 sy += p.y;
0755 sz += p.z;
0756 }
0757
0758 const double n = static_cast<double>(points.size());
0759 c.x = sx / n;
0760 c.y = sy / n;
0761 c.z = sz / n;
0762
0763 double sxx = 0.0;
0764 double syy = 0.0;
0765 double szz = 0.0;
0766 for (const Point& p : points)
0767 {
0768 const double dx = p.x - c.x;
0769 const double dy = p.y - c.y;
0770 const double dz = p.z - c.z;
0771 sxx += dx * dx;
0772 syy += dy * dy;
0773 szz += dz * dz;
0774 }
0775
0776 c.rms_x = std::sqrt(sxx / n);
0777 c.rms_y = std::sqrt(syy / n);
0778 c.rms_z = std::sqrt(szz / n);
0779 c.layer = points.front().layer;
0780 c.ok = std::isfinite(c.x) && std::isfinite(c.y) && std::isfinite(c.z);
0781 return c;
0782 }
0783
0784 int Tpc_PolyClusterizer::process_event(PHCompositeNode* topNode)
0785 {
0786 if (!m_assembledTracks || !m_clusters || !m_crossingDecisions || !m_garfield) { return Fun4AllReturnCodes::EVENT_OK;
0787 }
0788
0789 ActsGeometry* tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0790 if (!tGeometry)
0791 {
0792 std::cerr << Name() << "::process_event - missing ActsGeometry, using RMS fallback for errors" << std::endl;
0793 }
0794 m_clusters->Reset();
0795
0796 const unsigned int nassembled = m_assembledTracks->size();
0797 unsigned int nclusters = 0;
0798 std::map<TrkrDefs::hitsetkey, unsigned int> next_cluster_index_by_hitset;
0799
0800 for (int side = 0; side < 2; ++side)
0801 {
0802 for (unsigned int sector = 0; sector < 12; ++sector)
0803 {
0804 for (unsigned int iassembled = 0; iassembled < nassembled; ++iassembled)
0805 {
0806 const Tpc_AssembledTrack* assembled = m_assembledTracks->get_track(iassembled);
0807 if (!assembled) { continue;
0808 }
0809 if (assembled->get_side() != side) { continue;
0810 }
0811 if (assembled->get_first_sector() % 12U != sector) { continue;
0812 }
0813 const TpcCrossingDecision* crossing_decision = m_crossingDecisions->get_decision(assembled->get_track_id());
0814 if (!crossing_decision) { continue;
0815 }
0816 const unsigned char selected_tier = crossing_decision->get_selected_tier();
0817 if (selected_tier > m_maxAcceptedTier) { continue;
0818 }
0819 const short selected_crossing = crossing_decision->get_selected_crossing();
0820
0821
0822 std::map<TrkrDefs::hitsetkey, std::vector<Point>> points_by_hitset;
0823 for (unsigned int ih = 0; ih < assembled->size_hit_indices(); ++ih)
0824 {
0825 const Tpc_AssembledTrack::HitIndex hi = assembled->get_hit_index(ih);
0826 if (TpcDefs::getSide(hi.first) != static_cast<unsigned int>(side)) { continue;
0827 }
0828
0829 Point p;
0830 if (make_xyz_point(hi.first, hi.second, selected_crossing, p)) { points_by_hitset[p.hitsetkey].push_back(p);
0831 }
0832 }
0833 if (points_by_hitset.empty()) { continue;
0834 }
0835
0836 for (const auto& hitset_points : points_by_hitset)
0837 {
0838 const TrkrDefs::hitsetkey cluster_hitsetkey = hitset_points.first;
0839 const std::vector<Point>& points = hitset_points.second;
0840 const Centroid centroid = make_centroid(points);
0841 if (!centroid.ok) { continue;
0842 }
0843
0844 const unsigned int cluster_index = next_cluster_index_by_hitset[cluster_hitsetkey]++;
0845 const TrkrDefs::cluskey trkr_cluster_key = TrkrDefs::genClusKey(cluster_hitsetkey, cluster_index);
0846
0847 Tpc_PolyClusterv1* out = new Tpc_PolyClusterv1();
0848 out->set_event(m_event);
0849 out->set_cluster_id(m_clusters->size());
0850 out->set_source_assembled_track_id(assembled->get_track_id());
0851 out->set_trkr_cluster_key(trkr_cluster_key);
0852 out->set_side(side);
0853 out->set_centroid_x(centroid.x);
0854 out->set_centroid_y(centroid.y);
0855 out->set_centroid_z(centroid.z);
0856
0857 double phi_error = std::hypot(centroid.rms_x, centroid.rms_y);
0858 double z_error = std::fabs(centroid.rms_z);
0859 PHG4TpcGeom* layergeom = m_geomContainerTpc ? m_geomContainerTpc->GetLayerCellGeom(centroid.layer) : nullptr;
0860 if (layergeom && tGeometry)
0861 {
0862 double adc_sum = 0.0;
0863 double iphi_sum = 0.0;
0864 double iphi2_sum = 0.0;
0865 double t_sum = 0.0;
0866 double t2_sum = 0.0;
0867 int phibinhi = -1;
0868 int phibinlo = std::numeric_limits<int>::max();
0869 int tbinhi = -1;
0870 int tbinlo = std::numeric_limits<int>::max();
0871
0872 for (const Point& p : points)
0873 {
0874 if (p.adc <= 0.0) { continue;
0875 }
0876
0877 const int iphi = static_cast<int>(p.pad);
0878 const int it = static_cast<int>(p.tbin);
0879 const double adc = p.adc;
0880 phibinhi = std::max(iphi, phibinhi);
0881 phibinlo = std::min(iphi, phibinlo);
0882 tbinhi = std::max(it, tbinhi);
0883 tbinlo = std::min(it, tbinlo);
0884
0885 iphi_sum += static_cast<double>(iphi) * adc;
0886 iphi2_sum += square(static_cast<double>(iphi)) * adc;
0887
0888 const double t = layergeom->get_zcenter(it);
0889 t_sum += t * adc;
0890 t2_sum += square(t) * adc;
0891 adc_sum += adc;
0892 }
0893
0894 const double drift_velocity = tGeometry->get_drift_velocity();
0895 if (adc_sum > 0.0 && std::isfinite(drift_velocity))
0896 {
0897 const double radius = layergeom->get_radius();
0898 const double clusiphi = iphi_sum / adc_sum;
0899 const double clust = t_sum / adc_sum;
0900 const double phi_cov = std::max(0.0, (iphi2_sum / adc_sum - square(clusiphi)) * square(layergeom->get_phistep()));
0901 const double t_cov = std::max(0.0, t2_sum / adc_sum - square(clust));
0902 const double phi_err_square = (phibinhi == phibinlo) ?
0903 9.0 * (square(radius * layergeom->get_phistep()) / 12.0) :
0904 square(radius) * phi_cov / (adc_sum * 0.14);
0905 const double t_err_square = (tbinhi == tbinlo) ?
0906 9.0 * (square(layergeom->get_zstep()) / 12.0) :
0907 t_cov / (adc_sum * 0.14);
0908
0909 if (phi_err_square >= 0.0 && std::isfinite(phi_err_square)) { phi_error = std::sqrt(phi_err_square);
0910 }
0911 const double z_err_square = t_err_square * square(drift_velocity);
0912 if (z_err_square >= 0.0 && std::isfinite(z_err_square)) { z_error = std::sqrt(z_err_square);
0913 }
0914 }
0915 }
0916
0917 out->set_rms_x(phi_error);
0918 out->set_rms_y(0.0);
0919 out->set_rms_z(z_error);
0920
0921 const ClusterParameters params = make_cluster_parameters(points, centroid, static_cast<int>(points.front().side));
0922 out->set_adc(params.adc);
0923 out->set_phi_width(params.phi_width);
0924 out->set_time_width(params.time_width);
0925 out->set_phase(params.phase);
0926 for (const Point& p : points) { out->add_hit(p.hitsetkey, p.hitkey, p.x, p.y, p.z);
0927 }
0928 if (out->size_hits() == 0)
0929 {
0930 delete out;
0931 continue;
0932 }
0933 m_clusters->add_cluster(out);
0934 ++nclusters;
0935 }
0936 }
0937 }
0938 }
0939
0940 if (Verbosity() > 0)
0941 {
0942 std::cout << Name() << "::process_event - event " << m_event
0943 << " assembled_tracks=" << nassembled
0944 << " poly_clusters=" << m_clusters->size()
0945 << " layer_clusters=" << nclusters << std::endl;
0946 }
0947
0948 ++m_event;
0949 return Fun4AllReturnCodes::EVENT_OK;
0950 }
