File indexing completed on 2026-08-31 08:21:41
0001 #include "Tpc_AssembledTrackDisplay.h"
0002
0003 #include "tpctrackreco/IdealPadMap.h"
0004 #include "tpctrackreco/Tpc_AssembledTrack.h"
0005 #include "tpctrackreco/Tpc_AssembledTrackContainer.h"
0006 #include "tpctrackreco/Tpc_FittingTools.h"
0007
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/getClass.h>
0011
0012 #include <trackbase/TpcDefs.h>
0013 #include <trackbase/TrkrHit.h>
0014 #include <trackbase/TrkrHitSet.h>
0015 #include <trackbase/TrkrHitSetContainer.h>
0016
0017 #include <TCanvas.h>
0018 #include <TColor.h>
0019 #include <TDirectory.h>
0020 #include <TFile.h>
0021 #include <TH3D.h>
0022 #include <TMath.h>
0023 #include <TPolyLine3D.h>
0024
0025 #include <algorithm>
0026 #include <cmath>
0027 #include <format>
0028 #include <iostream>
0029 #include <limits>
0030 #include <set>
0031 #include <string>
0032 #include <vector>
0033
0034 namespace
0035 {
0036 double wrap_phi(double phi)
0037 {
0038 while (phi > TMath::Pi())
0039 {
0040 phi -= 2.0 * TMath::Pi();
0041 }
0042 while (phi <= -TMath::Pi())
0043 {
0044 phi += 2.0 * TMath::Pi();
0045 }
0046 return phi;
0047 }
0048
0049 double unwrap_phi_near(const double phi, const double reference)
0050 {
0051 double out = phi;
0052 while (out - reference > TMath::Pi())
0053 {
0054 out -= 2.0 * TMath::Pi();
0055 }
0056 while (out - reference < -TMath::Pi())
0057 {
0058 out += 2.0 * TMath::Pi();
0059 }
0060 return out;
0061 }
0062
0063 bool is_good_number(const double x)
0064 {
0065 return std::isfinite(x) && std::fabs(x) < 1.0e30;
0066 }
0067
0068 int track_color(const unsigned int itrk)
0069 {
0070 static const int colors[] = {
0071 kRed + 1, kBlue + 1, kGreen + 2, kMagenta + 1, kCyan + 2,
0072 kOrange + 7, kViolet + 1, kAzure + 1, kPink + 7, kTeal + 3};
0073 return colors[itrk % (sizeof(colors) / sizeof(colors[0]))];
0074 }
0075
0076 unsigned long long make_unique_hit_id(const TrkrDefs::hitsetkey hsk,
0077 const TrkrDefs::hitkey hk)
0078 {
0079 return (static_cast<unsigned long long>(hsk) << 32U) |
0080 static_cast<unsigned long long>(hk);
0081 }
0082
0083 void style_fit_line_3d(TPolyLine3D* line, const int color)
0084 {
0085 if (!line)
0086 {
0087 return;
0088 }
0089 line->SetLineColor(color);
0090 line->SetLineWidth(4);
0091 line->SetLineStyle(1);
0092 }
0093
0094 double sagitta_model_derivative(const double xrot,
0095 const double x0,
0096 const double invR)
0097 {
0098 const double dx = xrot - x0;
0099 const double dx2 = dx * dx;
0100 const double invR2 = invR * invR;
0101 const double invR3 = invR2 * invR;
0102 const double invR5 = invR3 * invR2;
0103 return -invR * dx - 0.5 * invR3 * dx2 * dx - 0.375 * invR5 * dx2 * dx2 * dx;
0104 }
0105
0106 double sagitta_phi_at_radius(const double radius,
0107 const Tpc_FittingTools::SagittaFit& fit)
0108 {
0109 const double c = std::cos(fit.theta);
0110 const double s = std::sin(fit.theta);
0111 double yy = std::tan(fit.theta) * radius;
0112
0113 for (unsigned int iter = 0; iter < 25; ++iter)
0114 {
0115 const double xrot = c * radius + s * yy;
0116 const double yrot = -s * radius + c * yy;
0117 const double f = Tpc_FittingTools::sagittaModel(xrot, fit.S, fit.x0, fit.invR);
0118 const double g = yrot - f;
0119 const double df = sagitta_model_derivative(xrot, fit.x0, fit.invR);
0120 const double dg = c - df * s;
0121 if (std::fabs(dg) < 1.0e-12)
0122 {
0123 break;
0124 }
0125 const double step = g / dg;
0126 yy -= step;
0127 if (std::fabs(step) < 1.0e-10)
0128 {
0129 break;
0130 }
0131 }
0132
0133 return fit.b + yy;
0134 }
0135
0136
0137 struct RadiusSort
0138 {
0139 const std::vector<Tpc_AssembledTrackDisplay::HitPoint>* pts;
0140 explicit RadiusSort(const std::vector<Tpc_AssembledTrackDisplay::HitPoint>* p)
0141 : pts(p)
0142 {
0143 }
0144 bool operator()(unsigned int a, unsigned int b) const
0145 {
0146 const Tpc_AssembledTrackDisplay::HitPoint& pa = (*pts)[a];
0147 const Tpc_AssembledTrackDisplay::HitPoint& pb = (*pts)[b];
0148 if (pa.radius != pb.radius)
0149 {
0150 return pa.radius < pb.radius;
0151 }
0152 return pa.tbin < pb.tbin;
0153 }
0154 };
0155
0156 struct FitResult
0157 {
0158 FitResult()
0159 : ok(false)
0160 , use_sagitta(false)
0161 , phi_slope(0.0)
0162 , phi_intercept(0.0)
0163 , tbin_slope(0.0)
0164 , tbin_intercept(0.0)
0165 , rmin(0.0)
0166 , rmax(0.0)
0167 {
0168 }
0169
0170 bool ok;
0171 bool use_sagitta;
0172 double phi_slope;
0173 double phi_intercept;
0174 double tbin_slope;
0175 double tbin_intercept;
0176 Tpc_FittingTools::SagittaFit phi_sagitta;
0177 double rmin;
0178 double rmax;
0179 };
0180
0181
0182
0183 FitResult fit_assembled_track_points(const std::vector<Tpc_AssembledTrackDisplay::HitPoint>& pts,
0184 const int fit_mode,
0185 const double weight_power,
0186 const double weight_floor_frac)
0187 {
0188 FitResult result;
0189 if (pts.size() < 2)
0190 {
0191 return result;
0192 }
0193
0194 std::vector<unsigned int> order;
0195 order.reserve(pts.size());
0196 for (unsigned int i = 0; i < pts.size(); ++i)
0197 {
0198 order.push_back(i);
0199 }
0200 std::sort(order.begin(), order.end(), RadiusSort(&pts));
0201
0202 double max_adc = 0.0;
0203 for (const auto& p : pts)
0204 {
0205 max_adc = std::max(max_adc, static_cast<double>(p.adc));
0206 }
0207 if (max_adc <= 0.0)
0208 {
0209 max_adc = 1.0;
0210 }
0211
0212 std::vector<Tpc_FittingTools::FitPoint> radius_phi_points;
0213 std::vector<Tpc_FittingTools::FitPoint> radius_tbin_points;
0214 radius_phi_points.reserve(pts.size());
0215 radius_tbin_points.reserve(pts.size());
0216
0217 bool first = true;
0218 double phi_reference = 0.0;
0219 result.rmin = std::numeric_limits<double>::max();
0220 result.rmax = -std::numeric_limits<double>::max();
0221
0222 for (unsigned int io : order)
0223 {
0224 const Tpc_AssembledTrackDisplay::HitPoint& p = pts[io];
0225 if (!is_good_number(p.radius) || !is_good_number(p.global_phi))
0226 {
0227 continue;
0228 }
0229
0230 double phi = p.global_phi;
0231 if (first)
0232 {
0233 phi_reference = phi;
0234 first = false;
0235 }
0236 else
0237 {
0238 phi = unwrap_phi_near(phi, phi_reference);
0239 phi_reference = phi;
0240 }
0241
0242 const double w = Tpc_FittingTools::adcWeight(static_cast<double>(p.adc), max_adc,
0243 weight_power, weight_floor_frac);
0244 radius_phi_points.emplace_back(p.radius, phi, w);
0245 radius_tbin_points.emplace_back(p.radius, static_cast<double>(p.tbin), w);
0246 result.rmin = std::min(result.rmin, p.radius);
0247 result.rmax = std::max(result.rmax, p.radius);
0248 }
0249
0250 if (radius_phi_points.size() < 2 || radius_tbin_points.size() < 2)
0251 {
0252 return result;
0253 }
0254
0255 const Tpc_FittingTools::LineFit phi_fit = Tpc_FittingTools::fitLine(radius_phi_points);
0256 const Tpc_FittingTools::LineFit tbin_fit = Tpc_FittingTools::fitLine(radius_tbin_points);
0257 if (!phi_fit.ok || !tbin_fit.ok)
0258 {
0259 return result;
0260 }
0261
0262 result.ok = true;
0263 result.phi_slope = phi_fit.slope;
0264 result.phi_intercept = phi_fit.intercept;
0265 result.tbin_slope = tbin_fit.slope;
0266 result.tbin_intercept = tbin_fit.intercept;
0267
0268 if (fit_mode == Tpc_FittingTools::FIT_SAGITTA && radius_phi_points.size() >= 3)
0269 {
0270 result.phi_sagitta = Tpc_FittingTools::fitSagitta(radius_phi_points);
0271 result.use_sagitta = result.phi_sagitta.ok;
0272 }
0273
0274 return result;
0275 }
0276
0277 double fit_direction(const FitResult& fit)
0278 {
0279 if (!fit.ok)
0280 {
0281 return 0.0;
0282 }
0283 return fit.use_sagitta ? fit.phi_sagitta.theta : std::atan(fit.phi_slope);
0284 }
0285
0286 double fit_theta(const FitResult& fit)
0287 {
0288 if (!fit.ok)
0289 {
0290 return 0.0;
0291 }
0292 return std::atan(fit.tbin_slope);
0293 }
0294
0295 double fit_curvature(const FitResult& fit)
0296 {
0297 if (!fit.ok)
0298 {
0299 return 0.0;
0300 }
0301 return fit.use_sagitta ? fit.phi_sagitta.invR : 0.0;
0302 }
0303
0304 void flush_fit_segment(std::vector<TPolyLine3D*>& lines,
0305 std::vector<double>& tbin_values,
0306 std::vector<double>& phi_values,
0307 std::vector<double>& radius_values,
0308 const int color,
0309 const bool draw_xy)
0310 {
0311 if (tbin_values.size() < 2)
0312 {
0313 tbin_values.clear();
0314 phi_values.clear();
0315 radius_values.clear();
0316 return;
0317 }
0318
0319 TPolyLine3D* line = new TPolyLine3D(static_cast<int>(tbin_values.size()));
0320 for (unsigned int i = 0; i < tbin_values.size(); ++i)
0321 {
0322 if (draw_xy)
0323 {
0324 const double x = radius_values[i] * std::cos(phi_values[i]);
0325 const double y = radius_values[i] * std::sin(phi_values[i]);
0326 line->SetPoint(static_cast<int>(i), tbin_values[i], x, y);
0327 }
0328 else
0329 {
0330 line->SetPoint(static_cast<int>(i), tbin_values[i], phi_values[i], radius_values[i]);
0331 }
0332 }
0333 style_fit_line_3d(line, color);
0334 lines.push_back(line);
0335
0336 tbin_values.clear();
0337 phi_values.clear();
0338 radius_values.clear();
0339 }
0340
0341 void add_fit_lines(std::vector<TPolyLine3D*>& lines,
0342 const FitResult& fit,
0343 const int color,
0344 const bool draw_xy)
0345 {
0346 if (!fit.ok || fit.rmax <= fit.rmin)
0347 {
0348 return;
0349 }
0350
0351 const int npts = 120;
0352 std::vector<double> seg_tbin;
0353 std::vector<double> seg_phi;
0354 std::vector<double> seg_radius;
0355 double previous_phi = 0.0;
0356 bool have_previous = false;
0357
0358 for (int i = 0; i < npts; ++i)
0359 {
0360 const double f = static_cast<double>(i) / static_cast<double>(npts - 1);
0361 const double radius = fit.rmin + f * (fit.rmax - fit.rmin);
0362 double phi = fit.use_sagitta ? sagitta_phi_at_radius(radius, fit.phi_sagitta)
0363 : fit.phi_slope * radius + fit.phi_intercept;
0364 phi = wrap_phi(phi);
0365 const double tbin = fit.tbin_slope * radius + fit.tbin_intercept;
0366
0367 if (!is_good_number(phi) || !is_good_number(tbin) || !is_good_number(radius))
0368 {
0369 continue;
0370 }
0371
0372 const bool crosses_phi_edge = have_previous && std::fabs(phi - previous_phi) > TMath::Pi();
0373 if (crosses_phi_edge)
0374 {
0375 flush_fit_segment(lines, seg_tbin, seg_phi, seg_radius, color, draw_xy);
0376 }
0377
0378 seg_tbin.push_back(tbin);
0379 seg_phi.push_back(phi);
0380 seg_radius.push_back(radius);
0381 previous_phi = phi;
0382 have_previous = true;
0383 }
0384
0385 flush_fit_segment(lines, seg_tbin, seg_phi, seg_radius, color, draw_xy);
0386 }
0387 }
0388
0389 Tpc_AssembledTrackDisplay::HitPoint::HitPoint()
0390 : ok(false)
0391 , hitsetkey(0)
0392 , hitkey(0)
0393 , layer(0)
0394 , region(0)
0395 , sector(0)
0396 , side(0)
0397 , pad(0)
0398 , tbin(0)
0399 , adc(0)
0400 , global_phi(0.0)
0401 , radius(0.0)
0402 {
0403 }
0404
0405 Tpc_AssembledTrackDisplay::Tpc_AssembledTrackDisplay(const std::string& name,
0406 const std::string& outfilename,
0407 const std::string& trackNodeName,
0408 const unsigned int maxEventDisplays)
0409 : SubsysReco(name)
0410 , m_outfilename(outfilename)
0411 , m_trackNodeName(trackNodeName)
0412 , m_maxEventDisplays(maxEventDisplays)
0413 , m_evt(0)
0414 , m_eventsSaved(0)
0415 , m_outfile(nullptr)
0416 , m_tracks(nullptr)
0417 , m_hits(nullptr)
0418 , m_idealPadMap(new IdealPadMap())
0419 , m_fitMode(Tpc_FittingTools::FIT_SAGITTA)
0420 , m_fitWeightPower(1.0)
0421 , m_fitWeightFloorFrac(0.05)
0422 , m_plotMinModules(0)
0423 , m_plotMaxModules(999999)
0424 , m_plotMinDirection(-1.0e30)
0425 , m_plotMaxDirection(1.0e30)
0426 , m_plotMinTheta(-1.0e30)
0427 , m_plotMaxTheta(1.0e30)
0428 , m_plotMinCurvature(-1.0e30)
0429 , m_plotMaxCurvature(1.0e30)
0430 {
0431 }
0432
0433 Tpc_AssembledTrackDisplay::~Tpc_AssembledTrackDisplay()
0434 {
0435 delete m_idealPadMap;
0436 m_idealPadMap = nullptr;
0437 }
0438
0439 int Tpc_AssembledTrackDisplay::Init(PHCompositeNode* )
0440 {
0441
0442 m_outfile = new TFile(m_outfilename.c_str(), "RECREATE");
0443 if (!m_outfile || m_outfile->IsZombie())
0444 {
0445 std::cerr << "Tpc_AssembledTrackDisplay::Init - cannot open output file "
0446 << m_outfilename << std::endl;
0447 return Fun4AllReturnCodes::ABORTRUN;
0448 }
0449
0450 if (!m_idealPadMap)
0451 {
0452 m_idealPadMap = new IdealPadMap();
0453 }
0454 if (!m_idealPadMap->is_loaded() && m_idealPadMap->load_from_cdb(Verbosity()) != 0)
0455 {
0456 std::cerr << "Tpc_AssembledTrackDisplay::Init - failed to load IdealPadMap from CDB" << std::endl;
0457 return Fun4AllReturnCodes::ABORTRUN;
0458 }
0459
0460 m_outfile->mkdir("events");
0461 std::cout << "Tpc_AssembledTrackDisplay::Init - writing up to "
0462 << m_maxEventDisplays << " events to "
0463 << m_outfilename << std::endl;
0464
0465 return Fun4AllReturnCodes::EVENT_OK;
0466 }
0467
0468 int Tpc_AssembledTrackDisplay::process_event(PHCompositeNode* topNode)
0469 {
0470 ++m_evt;
0471
0472 if (!get_nodes(topNode))
0473 {
0474 return Fun4AllReturnCodes::EVENT_OK;
0475 }
0476 if (!m_outfile || m_eventsSaved >= m_maxEventDisplays)
0477 {
0478 return Fun4AllReturnCodes::EVENT_OK;
0479 }
0480
0481 const unsigned int ntracks = m_tracks ? m_tracks->size() : 0;
0482
0483 TDirectory* eventsTop = m_outfile->GetDirectory("events");
0484 if (!eventsTop)
0485 {
0486 eventsTop = m_outfile->mkdir("events");
0487 }
0488 eventsTop->cd();
0489
0490 TDirectory* eventDir = eventsTop->mkdir(std::format("event_{:06}", m_evt).c_str());
0491 if (!eventDir)
0492 {
0493 std::cerr << "Tpc_AssembledTrackDisplay::process_event - failed to create event directory" << std::endl;
0494 return Fun4AllReturnCodes::EVENT_OK;
0495 }
0496 eventDir->cd();
0497
0498 TH3D* h3[2] = {nullptr, nullptr};
0499 TH3D* h3xy[2] = {nullptr, nullptr};
0500 TH3D* h3_single[2] = {nullptr, nullptr};
0501 TH3D* h3xy_single[2] = {nullptr, nullptr};
0502
0503 std::vector<TPolyLine3D*> fit_lines_tpr[2];
0504 std::vector<TPolyLine3D*> fit_lines_txy[2];
0505 std::vector<TPolyLine3D*> fit_lines_tpr_single[2];
0506 std::vector<TPolyLine3D*> fit_lines_txy_single[2];
0507
0508 std::set<unsigned long long> filled_hit_ids[2];
0509 std::set<unsigned long long> filled_hit_ids_single[2];
0510
0511 for (unsigned int side = 0; side < 2; ++side)
0512 {
0513 h3[side] = new TH3D(std::format("h3_evt{:06}_tpc_assembledtrack_hits_side{}", m_evt, side).c_str(),
0514 std::format("event {} side {} assembled tracks;timebin;global #phi;radius [cm]", m_evt, side).c_str(),
0515 512, -0.5, 511.5,
0516 720, -TMath::Pi(), TMath::Pi(),
0517 100, 30, 80.0);
0518 h3[side]->SetStats(false);
0519 h3[side]->SetDirectory(nullptr);
0520
0521 h3xy[side] = new TH3D(std::format("h3_evt{:06}_tpc_assembledtrack_hits_xy_side{}", m_evt, side).c_str(),
0522 std::format("event {} side {} assembled tracks;timebin;x [cm];y [cm]", m_evt, side).c_str(),
0523 512, -0.5, 511.5,
0524 160, -80.0, 80.0,
0525 160, -80.0, 80.0);
0526 h3xy[side]->SetStats(false);
0527 h3xy[side]->SetDirectory(nullptr);
0528
0529 h3_single[side] = new TH3D(std::format("h3_evt{:06}_single_module_tpc_assembledtrack_hits_side{}", m_evt, side).c_str(),
0530 std::format("event {} side {} single-module assembled tracks;timebin;global #phi;radius [cm]", m_evt, side).c_str(),
0531 512, -0.5, 511.5,
0532 720, -TMath::Pi(), TMath::Pi(),
0533 100, 30, 80.0);
0534 h3_single[side]->SetStats(false);
0535 h3_single[side]->SetDirectory(nullptr);
0536
0537 h3xy_single[side] = new TH3D(std::format("h3_evt{:06}_single_module_tpc_assembledtrack_hits_xy_side{}", m_evt, side).c_str(),
0538 std::format("event {} side {} single-module assembled tracks;timebin;x [cm];y [cm]", m_evt, side).c_str(),
0539 512, -0.5, 511.5,
0540 160, -80.0, 80.0,
0541 160, -80.0, 80.0);
0542 h3xy_single[side]->SetStats(false);
0543 h3xy_single[side]->SetDirectory(nullptr);
0544 }
0545
0546 for (unsigned int itrk = 0; itrk < ntracks; ++itrk)
0547 {
0548 const Tpc_AssembledTrack* trk = m_tracks->get_track(itrk);
0549 if (!trk)
0550 {
0551 continue;
0552 }
0553
0554 const int trk_side = trk->get_side();
0555 if (trk_side < 0 || trk_side > 1)
0556 {
0557 continue;
0558 }
0559 const unsigned int side = static_cast<unsigned int>(trk_side);
0560 const bool single_module_track = (trk->get_nsegments() == 1U);
0561
0562 const unsigned int nmodules = trk->get_nsegments();
0563 if (nmodules < m_plotMinModules || nmodules > m_plotMaxModules)
0564 {
0565 continue;
0566 }
0567
0568 std::vector<HitPoint> pts;
0569 pts.reserve(trk->size_hit_indices());
0570
0571 for (unsigned int ih = 0; ih < trk->size_hit_indices(); ++ih)
0572 {
0573 const Tpc_AssembledTrack::HitIndex idx = trk->get_hit_index(ih);
0574 const HitPoint p = make_hit_point(idx.first, idx.second);
0575 if (!p.ok)
0576 {
0577 continue;
0578 }
0579 pts.push_back(p);
0580 }
0581
0582 const FitResult fit = fit_assembled_track_points(pts, m_fitMode, m_fitWeightPower, m_fitWeightFloorFrac);
0583 if (!fit.ok)
0584 {
0585 continue;
0586 }
0587
0588 const double direction = fit_direction(fit);
0589 const double theta = fit_theta(fit);
0590 const double curvature = fit_curvature(fit);
0591 if (direction < m_plotMinDirection || direction > m_plotMaxDirection)
0592 {
0593 continue;
0594 }
0595 if (theta < m_plotMinTheta || theta > m_plotMaxTheta)
0596 {
0597 continue;
0598 }
0599 if (curvature < m_plotMinCurvature || curvature > m_plotMaxCurvature)
0600 {
0601 continue;
0602 }
0603
0604 for (const auto& p : pts)
0605 {
0606 const unsigned long long uid = make_unique_hit_id(p.hitsetkey, p.hitkey);
0607 const double phi = wrap_phi(p.global_phi);
0608 const double x = p.radius * std::cos(p.global_phi);
0609 const double y = p.radius * std::sin(p.global_phi);
0610
0611 if (filled_hit_ids[side].insert(uid).second)
0612 {
0613 h3[side]->Fill(static_cast<double>(p.tbin), phi, p.radius, static_cast<double>(p.adc));
0614 h3xy[side]->Fill(static_cast<double>(p.tbin), x, y, static_cast<double>(p.adc));
0615 }
0616
0617 if (single_module_track && filled_hit_ids_single[side].insert(uid).second)
0618 {
0619 h3_single[side]->Fill(static_cast<double>(p.tbin), phi, p.radius, static_cast<double>(p.adc));
0620 h3xy_single[side]->Fill(static_cast<double>(p.tbin), x, y, static_cast<double>(p.adc));
0621 }
0622 }
0623
0624 const int color = track_color(itrk);
0625 add_fit_lines(fit_lines_tpr[side], fit, color, false);
0626 add_fit_lines(fit_lines_txy[side], fit, color, true);
0627
0628 if (single_module_track)
0629 {
0630 add_fit_lines(fit_lines_tpr_single[side], fit, color, false);
0631 add_fit_lines(fit_lines_txy_single[side], fit, color, true);
0632 }
0633 }
0634
0635 for (unsigned int side = 0; side < 2; ++side)
0636 {
0637 TCanvas* c3 = new TCanvas(std::format("c3_evt{:06}_timebin_phi_radius_fits_side{}", m_evt, side).c_str(),
0638 std::format("event {} side {} assembled-track hits and display fits", m_evt, side).c_str(),
0639 1200, 900);
0640 h3[side]->Draw("BOX2Z");
0641 for (auto& iline : fit_lines_tpr[side])
0642 {
0643 if (iline)
0644 {
0645 iline->Draw("same");
0646 }
0647 }
0648 c3->Modified();
0649 c3->Update();
0650 c3->Write();
0651
0652 TCanvas* c3xy = new TCanvas(std::format("c3_evt{:06}_timebin_x_y_fits_side{}", m_evt, side).c_str(),
0653 std::format("event {} side {} assembled-track hits and display fits", m_evt, side).c_str(),
0654 1200, 900);
0655 h3xy[side]->Draw("BOX2Z");
0656 for (auto& iline : fit_lines_txy[side])
0657 {
0658 if (iline)
0659 {
0660 iline->Draw("same");
0661 }
0662 }
0663 c3xy->Modified();
0664 c3xy->Update();
0665 c3xy->Write();
0666
0667 TCanvas* c3_single = new TCanvas(std::format("c3_evt{:06}_timebin_phi_radius_fits_single_module_side{}", m_evt, side).c_str(),
0668 std::format("event {} side {} single-module assembled-track hits and display fits", m_evt, side).c_str(),
0669 1200, 900);
0670 h3_single[side]->Draw("BOX2Z");
0671 for (auto& iline : fit_lines_tpr_single[side])
0672 {
0673 if (iline)
0674 {
0675 iline->Draw("same");
0676 }
0677 }
0678 c3_single->Modified();
0679 c3_single->Update();
0680 c3_single->Write();
0681
0682
0683
0684
0685
0686
0687
0688
0689
0690
0691
0692
0693
0694 }
0695
0696 std::cout << "Tpc_AssembledTrackDisplay - saved event " << m_evt
0697 << " with " << ntracks << " assembled tracks" << std::endl;
0698
0699 ++m_eventsSaved;
0700 return Fun4AllReturnCodes::EVENT_OK;
0701 }
0702
0703 int Tpc_AssembledTrackDisplay::End(PHCompositeNode* )
0704 {
0705 if (m_outfile)
0706 {
0707 m_outfile->Close();
0708 delete m_outfile;
0709 m_outfile = nullptr;
0710 }
0711
0712 std::cout << "Tpc_AssembledTrackDisplay::End - events seen: " << m_evt
0713 << ", events written: " << m_eventsSaved
0714 << ", output file: " << m_outfilename << std::endl;
0715
0716 return Fun4AllReturnCodes::EVENT_OK;
0717 }
0718
0719 bool Tpc_AssembledTrackDisplay::get_nodes(PHCompositeNode* topNode)
0720 {
0721 m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0722 m_tracks = findNode::getClass<Tpc_AssembledTrackContainer>(topNode, m_trackNodeName);
0723
0724 if (!m_tracks)
0725 {
0726 const char* candidate_names[] = {
0727 "TPC_ASSEMBLEDTRACKS",
0728 "Tpc_AssembledTracks",
0729 "Tpc_AssembledTrackContainer",
0730 "TPC_ASSEMBLEDTRACKCONTAINER",
0731 "TPC_ASSEMBLEDTRACKS_CONTAINER"};
0732
0733 for (unsigned int i = 0;
0734 i < sizeof(candidate_names) / sizeof(candidate_names[0]) && !m_tracks;
0735 ++i)
0736 {
0737 m_tracks = findNode::getClass<Tpc_AssembledTrackContainer>(topNode, candidate_names[i]);
0738 if (m_tracks)
0739 {
0740 m_trackNodeName = candidate_names[i];
0741 }
0742 }
0743 }
0744
0745 if (!m_idealPadMap)
0746 {
0747 m_idealPadMap = new IdealPadMap();
0748 }
0749
0750 if (!m_tracks)
0751 {
0752 std::cerr << "Tpc_AssembledTrackDisplay - could not find Tpc_AssembledTrackContainer node" << std::endl;
0753 return false;
0754 }
0755
0756 if (!m_hits)
0757 {
0758 std::cerr << "Tpc_AssembledTrackDisplay - missing TRKR_HITSET" << std::endl;
0759 return false;
0760 }
0761
0762 if (!m_idealPadMap || !m_idealPadMap->is_loaded())
0763 {
0764 std::cerr << "Tpc_AssembledTrackDisplay - IdealPadMap is not loaded" << std::endl;
0765 return false;
0766 }
0767
0768 return true;
0769 }
0770
0771 Tpc_AssembledTrackDisplay::HitPoint
0772 Tpc_AssembledTrackDisplay::make_hit_point(const TrkrDefs::hitsetkey hsk,
0773 const TrkrDefs::hitkey hk) const
0774 {
0775 HitPoint p;
0776
0777 TrkrHitSet* hitset = m_hits ? m_hits->findHitSet(hsk) : nullptr;
0778 if (!hitset)
0779 {
0780 return p;
0781 }
0782
0783 TrkrHit* hit = hitset->getHit(hk);
0784 if (!hit)
0785 {
0786 return p;
0787 }
0788
0789 p.hitsetkey = hsk;
0790 p.hitkey = hk;
0791 p.layer = TrkrDefs::getLayer(hsk);
0792 p.side = static_cast<int>(TpcDefs::getSide(hsk));
0793 p.pad = TpcDefs::getPad(hk);
0794 p.tbin = TpcDefs::getTBin(hk);
0795 p.adc = hit->getAdc();
0796
0797 if (p.layer < 7 || p.layer > 54)
0798 {
0799 return p;
0800 }
0801 p.region = static_cast<unsigned int>((p.layer - 7) / 16);
0802
0803 if (!m_idealPadMap)
0804 {
0805 return p;
0806 }
0807 const unsigned int pads_per_sector = m_idealPadMap->get_pads_per_sector_for_layer(p.layer);
0808 if (pads_per_sector == 0U)
0809 {
0810 return p;
0811 }
0812
0813 p.sector = (p.pad / pads_per_sector) % 12U;
0814 p.radius = m_idealPadMap->get_radius(p.layer);
0815 p.global_phi = wrap_phi(m_idealPadMap->get_phi(static_cast<unsigned int>(p.side), p.layer, p.pad));
0816
0817 if (!is_good_number(p.radius) || !is_good_number(p.global_phi))
0818 {
0819 return p;
0820 }
0821
0822 p.ok = true;
0823 return p;
0824 }
