File indexing completed on 2026-08-31 08:21:41
0001 #include "Tpc_ModuleTrackDisplay.h"
0002
0003 #include "tpctrackreco/IdealPadMap.h"
0004 #include "tpctrackreco/Tpc_FittingTools.h"
0005 #include "tpctrackreco/Tpc_ModuleTrack.h"
0006 #include "tpctrackreco/Tpc_ModuleTrackContainer.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 <TDirectory.h>
0019 #include <TFile.h>
0020 #include <TGraph.h>
0021 #include <TH3D.h>
0022 #include <TMath.h>
0023 #include <TMultiGraph.h>
0024 #include <TPolyLine3D.h>
0025
0026 #include <algorithm>
0027 #include <cmath>
0028 #include <format>
0029 #include <iostream>
0030 #include <map>
0031 #include <set>
0032 #include <string>
0033 #include <utility>
0034 #include <vector>
0035
0036 namespace
0037 {
0038 int track_color(const unsigned int itrk)
0039 {
0040 static const int colors[] = {
0041 kRed + 1, kBlue + 1, kGreen + 2, kMagenta + 1, kCyan + 2,
0042 kOrange + 7, kViolet + 1, kAzure + 1, kPink + 7, kTeal + 3};
0043 return colors[itrk % (sizeof(colors) / sizeof(colors[0]))];
0044 }
0045
0046 int layer_to_module(const unsigned int layer)
0047 {
0048 if (layer < 7 || layer > 54)
0049 {
0050 return -1;
0051 }
0052 const int module = static_cast<int>((layer - 7) / 16);
0053 return (module >= 0 && module < 3) ? module : -1;
0054 }
0055
0056 unsigned long long make_unique_hit_id(const TrkrDefs::hitsetkey hsk,
0057 const TrkrDefs::hitkey hk)
0058 {
0059 return (static_cast<unsigned long long>(hsk) << 32U) |
0060 static_cast<unsigned long long>(hk);
0061 }
0062
0063 bool is_good_number(const double x)
0064 {
0065 return std::isfinite(x) && std::fabs(x) < 1.0e30;
0066 }
0067
0068 double unwrap_phi_near(const double phi, const double reference)
0069 {
0070 double out = phi;
0071 while (out - reference > TMath::Pi())
0072 {
0073 out -= 2.0 * TMath::Pi();
0074 }
0075 while (out - reference < -TMath::Pi())
0076 {
0077 out += 2.0 * TMath::Pi();
0078 }
0079 return out;
0080 }
0081
0082 bool sagitta_y_at_x(const double x,
0083 const double mean_y,
0084 const Tpc_FittingTools::SagittaFit& fit,
0085 double& y_out)
0086 {
0087 if (!fit.ok || !is_good_number(x) || !is_good_number(mean_y))
0088 {
0089 return false;
0090 }
0091 if (!is_good_number(fit.S) || !is_good_number(fit.x0) ||
0092 !is_good_number(fit.invR) || !is_good_number(fit.theta) ||
0093 !is_good_number(fit.b))
0094 {
0095 return false;
0096 }
0097
0098 const double c = std::cos(fit.theta);
0099 const double st = std::sin(fit.theta);
0100 const double ymin = mean_y - 0.8;
0101 const double ymax = mean_y + 0.8;
0102 const int nscan = 240;
0103
0104 bool have_best = false;
0105 double best_y = 0.0;
0106 double best_dist = 1.0e99;
0107
0108 double yprev = ymin;
0109 double xprev_rot = c * x + st * (yprev - fit.b);
0110 double yprev_rot = -st * x + c * (yprev - fit.b);
0111 double gprev = yprev_rot - Tpc_FittingTools::sagittaModel(xprev_rot, fit.S, fit.x0, fit.invR);
0112
0113 double best_abs_g = std::fabs(gprev);
0114 double best_abs_y = yprev;
0115
0116 for (int iscan = 1; iscan <= nscan; ++iscan)
0117 {
0118 const double ycur = ymin + (ymax - ymin) * static_cast<double>(iscan) / static_cast<double>(nscan);
0119 const double xcur_rot = c * x + st * (ycur - fit.b);
0120 const double ycur_rot = -st * x + c * (ycur - fit.b);
0121 const double gcur = ycur_rot - Tpc_FittingTools::sagittaModel(xcur_rot, fit.S, fit.x0, fit.invR);
0122
0123 if (std::fabs(gcur) < best_abs_g)
0124 {
0125 best_abs_g = std::fabs(gcur);
0126 best_abs_y = ycur;
0127 }
0128
0129 const bool bracket =
0130 gprev == 0.0 || gcur == 0.0 ||
0131 (gprev < 0.0 && gcur > 0.0) ||
0132 (gprev > 0.0 && gcur < 0.0);
0133
0134 if (bracket)
0135 {
0136 double ya = yprev;
0137 double yb = ycur;
0138 double ga = gprev;
0139
0140 for (int ib = 0; ib < 50; ++ib)
0141 {
0142 const double ym = 0.5 * (ya + yb);
0143 const double xm_rot = c * x + st * (ym - fit.b);
0144 const double ym_rot = -st * x + c * (ym - fit.b);
0145 const double gm = ym_rot - Tpc_FittingTools::sagittaModel(xm_rot, fit.S, fit.x0, fit.invR);
0146
0147 if ((ga < 0.0 && gm <= 0.0) || (ga > 0.0 && gm >= 0.0))
0148 {
0149 ya = ym;
0150 ga = gm;
0151 }
0152 else
0153 {
0154 yb = ym;
0155 }
0156 }
0157
0158 const double y = 0.5 * (ya + yb);
0159 const double dist = std::fabs(y - mean_y);
0160 if (!have_best || dist < best_dist)
0161 {
0162 have_best = true;
0163 best_dist = dist;
0164 best_y = y;
0165 }
0166 }
0167
0168 yprev = ycur;
0169 gprev = gcur;
0170 }
0171
0172 if (have_best)
0173 {
0174 y_out = best_y;
0175 return true;
0176 }
0177
0178 if (best_abs_g < 2.0e-2)
0179 {
0180 y_out = best_abs_y;
0181 return true;
0182 }
0183
0184 return false;
0185 }
0186
0187 void style_hit_graph(TGraph* g, const int color)
0188 {
0189 if (!g)
0190 {
0191 return;
0192 }
0193 g->SetMarkerStyle(20);
0194 g->SetMarkerSize(0.8);
0195 g->SetMarkerColor(color);
0196 g->SetLineColor(color);
0197 g->SetLineWidth(2);
0198 }
0199
0200 void style_fit_graph(TGraph* g, const int color)
0201 {
0202 if (!g)
0203 {
0204 return;
0205 }
0206 g->SetMarkerStyle(1);
0207 g->SetMarkerSize(0.0);
0208 g->SetLineColor(color);
0209 g->SetLineWidth(3);
0210 g->SetLineStyle(2);
0211 }
0212
0213 void style_fit_line_3d(TPolyLine3D* line, const int color)
0214 {
0215 if (!line)
0216 {
0217 return;
0218 }
0219 line->SetLineColor(color);
0220 line->SetLineWidth(4);
0221 line->SetLineStyle(1);
0222 }
0223
0224
0225 struct GroupKey
0226 {
0227 GroupKey()
0228 : side(-1)
0229 , sector(-1)
0230 , module(-1)
0231 {
0232 }
0233 GroupKey(const int s, const int sec, const int mod)
0234 : side(s)
0235 , sector(sec)
0236 , module(mod)
0237 {
0238 }
0239
0240 int side;
0241 int sector;
0242 int module;
0243
0244 bool operator<(const GroupKey& rhs) const
0245 {
0246 if (side != rhs.side)
0247 {
0248 return side < rhs.side;
0249 }
0250 if (sector != rhs.sector)
0251 {
0252 return sector < rhs.sector;
0253 }
0254 return module < rhs.module;
0255 }
0256 };
0257
0258 struct FitResult
0259 {
0260 FitResult()
0261 : ok(false)
0262 , use_sagitta(false)
0263 , phi_slope(0.0)
0264 , phi_intercept(0.0)
0265 , tbin_slope(0.0)
0266 , tbin_intercept(0.0)
0267 , mean_phi(0.0)
0268 , rmin(0.0)
0269 , rmax(0.0)
0270 {
0271 }
0272
0273 bool ok;
0274 bool use_sagitta;
0275 double phi_slope;
0276 double phi_intercept;
0277 double tbin_slope;
0278 double tbin_intercept;
0279 Tpc_FittingTools::SagittaFit phi_sagitta;
0280 double mean_phi;
0281 double rmin;
0282 double rmax;
0283 };
0284
0285 struct HardwareFitResult
0286 {
0287 HardwareFitResult()
0288 : ok(false)
0289 , use_sagitta(false)
0290 , pad_slope(0.0)
0291 , pad_intercept(0.0)
0292 , tbin_slope(0.0)
0293 , tbin_intercept(0.0)
0294 , mean_pad(0.0)
0295 , layer_min(0.0)
0296 , layer_max(0.0)
0297 {
0298 }
0299
0300 bool ok;
0301 bool use_sagitta;
0302 double pad_slope;
0303 double pad_intercept;
0304 double tbin_slope;
0305 double tbin_intercept;
0306 Tpc_FittingTools::SagittaFit pad_sagitta;
0307 double mean_pad;
0308 double layer_min;
0309 double layer_max;
0310 };
0311
0312 struct GraphBundle
0313 {
0314 GraphBundle()
0315 : mg_phi_radius_hits(nullptr)
0316 , mg_tbin_radius_hits(nullptr)
0317 , mg_tbin_phi_hits(nullptr)
0318 , mg_phi_radius_fits(nullptr)
0319 , mg_tbin_radius_fits(nullptr)
0320 , mg_tbin_phi_fits(nullptr)
0321 , h3_hardware_hits(nullptr)
0322 , h3_adc_hits(nullptr)
0323 , h3_adc_unassociated_hits(nullptr)
0324 , c3_hardware_hits_fits(nullptr)
0325 , c3_adc_hits_fits(nullptr)
0326 , c3_adc_unassociated_hits(nullptr)
0327 , phi_reference(0.0)
0328 {
0329 }
0330
0331 TMultiGraph* mg_phi_radius_hits;
0332 TMultiGraph* mg_tbin_radius_hits;
0333 TMultiGraph* mg_tbin_phi_hits;
0334
0335 TMultiGraph* mg_phi_radius_fits;
0336 TMultiGraph* mg_tbin_radius_fits;
0337 TMultiGraph* mg_tbin_phi_fits;
0338
0339
0340 TH3D* h3_hardware_hits;
0341
0342
0343 TH3D* h3_adc_hits;
0344 TH3D* h3_adc_unassociated_hits;
0345
0346 TCanvas* c3_hardware_hits_fits;
0347 TCanvas* c3_adc_hits_fits;
0348 TCanvas* c3_adc_unassociated_hits;
0349
0350 double phi_reference;
0351
0352 std::vector<TPolyLine3D*> hardware_fit_lines_3d;
0353 std::vector<std::string> hardware_fit_line_names_3d;
0354
0355 std::vector<TPolyLine3D*> fit_lines_3d;
0356 std::vector<std::string> fit_line_names_3d;
0357
0358 std::set<unsigned long long> filled_hit_ids;
0359 };
0360
0361
0362
0363 void make_bundle_graphs(GraphBundle& b,
0364 const GroupKey& key,
0365 const unsigned int evt,
0366 const IdealPadMap* idealPadMap)
0367 {
0368 if (!idealPadMap)
0369 {
0370 return;
0371 }
0372
0373 const std::string tag = std::format("s{}_sec{:02}_mod{}", key.side, key.sector, key.module);
0374
0375 b.mg_phi_radius_hits = new TMultiGraph();
0376 b.mg_phi_radius_hits->SetName(std::format("mg_{}_phi_radius_hits", tag).c_str());
0377 b.mg_phi_radius_hits->SetTitle(std::format("event {} side {} sector {} module {} hits;#phi;radius [cm]",
0378 evt, key.side, key.sector, key.module)
0379 .c_str());
0380
0381 b.mg_tbin_radius_hits = new TMultiGraph();
0382 b.mg_tbin_radius_hits->SetName(std::format("mg_{}_tbin_radius_hits", tag).c_str());
0383 b.mg_tbin_radius_hits->SetTitle(std::format("event {} side {} sector {} module {} hits;timebin;radius [cm]",
0384 evt, key.side, key.sector, key.module)
0385 .c_str());
0386
0387 b.mg_tbin_phi_hits = new TMultiGraph();
0388 b.mg_tbin_phi_hits->SetName(std::format("mg_{}_tbin_phi_hits", tag).c_str());
0389 b.mg_tbin_phi_hits->SetTitle(std::format("event {} side {} sector {} module {} hits;timebin;#phi",
0390 evt, key.side, key.sector, key.module)
0391 .c_str());
0392
0393 b.mg_phi_radius_fits = new TMultiGraph();
0394 b.mg_phi_radius_fits->SetName(std::format("mg_{}_phi_radius_fits", tag).c_str());
0395 b.mg_phi_radius_fits->SetTitle(std::format("event {} side {} sector {} module {} display fits;#phi;radius [cm]",
0396 evt, key.side, key.sector, key.module)
0397 .c_str());
0398
0399 b.mg_tbin_radius_fits = new TMultiGraph();
0400 b.mg_tbin_radius_fits->SetName(std::format("mg_{}_tbin_radius_fits", tag).c_str());
0401 b.mg_tbin_radius_fits->SetTitle(std::format("event {} side {} sector {} module {} display fits;timebin;radius [cm]",
0402 evt, key.side, key.sector, key.module)
0403 .c_str());
0404
0405 b.mg_tbin_phi_fits = new TMultiGraph();
0406 b.mg_tbin_phi_fits->SetName(std::format("mg_{}_tbin_phi_fits", tag).c_str());
0407 b.mg_tbin_phi_fits->SetTitle(std::format("event {} side {} sector {} module {} display fits;timebin;#phi",
0408 evt, key.side, key.sector, key.module)
0409 .c_str());
0410
0411 const int layer_min = 7 + 16 * key.module;
0412 const int layer_max = layer_min + 15;
0413 const double radius_min = idealPadMap->get_radius(static_cast<unsigned int>(layer_min));
0414 const double radius_max = idealPadMap->get_radius(static_cast<unsigned int>(layer_max));
0415 const double radius_min_guess = std::min(radius_min, radius_max) - 1.0;
0416 const double radius_max_guess = std::max(radius_min, radius_max) + 1.0;
0417
0418 const unsigned int pads_per_sector =
0419 idealPadMap->get_pads_per_sector(static_cast<unsigned int>(key.module));
0420
0421 if (pads_per_sector < 2U)
0422 {
0423 return;
0424 }
0425
0426 const unsigned int nPads = pads_per_sector;
0427 const int pad_min = static_cast<int>(nPads) * key.sector;
0428 const int pad_max = pad_min + static_cast<int>(nPads) - 1;
0429 const unsigned int layer_ref = static_cast<unsigned int>(layer_min);
0430
0431
0432 const double phi_first =
0433 idealPadMap->get_phi(key.side,
0434 static_cast<unsigned int>(key.sector),
0435 layer_ref,
0436 0);
0437
0438 const double phi_last_wrapped =
0439 idealPadMap->get_phi(key.side,
0440 static_cast<unsigned int>(key.sector),
0441 layer_ref,
0442 pads_per_sector - 1);
0443
0444 const double phi_last = unwrap_phi_near(phi_last_wrapped, phi_first);
0445 const double dphi = (phi_last - phi_first) / static_cast<double>(pads_per_sector - 1);
0446
0447 double phi_min = phi_first - 0.5 * dphi;
0448 double phi_max = phi_last + 0.5 * dphi;
0449 if (phi_max < phi_min)
0450 {
0451 std::swap(phi_min, phi_max);
0452 }
0453 b.phi_reference = 0.5 * (phi_min + phi_max);
0454
0455 b.h3_hardware_hits = new TH3D(std::format("h3_{}_hardware_hits", tag).c_str(),
0456 std::format("event {} side {} sector {} module {} associated hits;timebin;pad;layer",
0457 evt, key.side, key.sector, key.module)
0458 .c_str(),
0459 512, -0.5, 511.5,
0460 static_cast<int>(nPads), static_cast<double>(pad_min) - 0.5, static_cast<double>(pad_max) + 0.5,
0461 16, static_cast<double>(layer_min) - 0.5, static_cast<double>(layer_max) + 0.5);
0462 b.h3_hardware_hits->SetStats(false);
0463
0464 b.h3_adc_hits = new TH3D(std::format("h3_{}_adc_hits", tag).c_str(),
0465 std::format("event {} side {} sector {} module {} associated hits;timebin;#phi;radius [cm]",
0466 evt, key.side, key.sector, key.module)
0467 .c_str(),
0468 512, -0.5, 511.5,
0469 static_cast<int>(nPads), phi_min, phi_max,
0470 20, radius_min_guess, radius_max_guess);
0471 b.h3_adc_hits->SetStats(false);
0472
0473 b.h3_adc_unassociated_hits = new TH3D(std::format("h3_{}_adc_unassociated_hits", tag).c_str(),
0474 std::format("event {} side {} sector {} module {} unassociated hits;timebin;#phi;radius [cm]",
0475 evt, key.side, key.sector, key.module)
0476 .c_str(),
0477 512, -0.5, 511.5,
0478 static_cast<int>(nPads), phi_min, phi_max,
0479 16, radius_min_guess, radius_max_guess);
0480 b.h3_adc_unassociated_hits->SetStats(false);
0481 }
0482
0483 GraphBundle& get_bundle(std::map<GroupKey, GraphBundle>& bundles,
0484 const GroupKey& key,
0485 const unsigned int evt,
0486 const IdealPadMap* idealPadMap)
0487 {
0488 std::map<GroupKey, GraphBundle>::iterator it = bundles.find(key);
0489 if (it == bundles.end())
0490 {
0491 const std::pair<std::map<GroupKey, GraphBundle>::iterator, bool> inserted =
0492 bundles.insert(std::make_pair(key, GraphBundle()));
0493 it = inserted.first;
0494 make_bundle_graphs(it->second, key, evt, idealPadMap);
0495 }
0496 return it->second;
0497 }
0498
0499 HardwareFitResult fit_hardware_track_points(const std::vector<Tpc_ModuleTrackDisplay::HitPoint>& pts,
0500 const int fit_mode,
0501 const double weight_power,
0502 const double weight_floor_frac)
0503 {
0504 HardwareFitResult result;
0505 if (pts.size() < 2)
0506 {
0507 return result;
0508 }
0509
0510 double max_adc = 0.0;
0511 for (const auto& p : pts)
0512 {
0513 max_adc = std::max(max_adc, static_cast<double>(p.adc));
0514 }
0515 if (max_adc <= 0.0)
0516 {
0517 max_adc = 1.0;
0518 }
0519
0520 std::vector<Tpc_FittingTools::FitPoint> layer_pad_points;
0521 std::vector<Tpc_FittingTools::FitPoint> layer_tbin_points;
0522 layer_pad_points.reserve(pts.size());
0523 layer_tbin_points.reserve(pts.size());
0524
0525 result.layer_min = static_cast<double>(pts.front().layer);
0526 result.layer_max = static_cast<double>(pts.front().layer);
0527 result.mean_pad = 0.0;
0528
0529 for (const auto& p : pts)
0530 {
0531 const double layer = static_cast<double>(p.layer);
0532 const double pad = static_cast<double>(p.pad);
0533 const double tbin = static_cast<double>(p.tbin);
0534 const double w = Tpc_FittingTools::adcWeight(static_cast<double>(p.adc), max_adc,
0535 weight_power, weight_floor_frac);
0536
0537 layer_pad_points.emplace_back(layer, pad, w);
0538 layer_tbin_points.emplace_back(layer, tbin, w);
0539 result.mean_pad += pad;
0540 result.layer_min = std::min(result.layer_min, layer);
0541 result.layer_max = std::max(result.layer_max, layer);
0542 }
0543
0544 result.mean_pad /= static_cast<double>(layer_pad_points.size());
0545
0546 const Tpc_FittingTools::LineFit pad_fit = Tpc_FittingTools::fitLine(layer_pad_points);
0547 const Tpc_FittingTools::LineFit tbin_fit = Tpc_FittingTools::fitLine(layer_tbin_points);
0548
0549 if (!pad_fit.ok || !tbin_fit.ok)
0550 {
0551 return result;
0552 }
0553
0554 result.ok = true;
0555 result.pad_slope = pad_fit.slope;
0556 result.pad_intercept = pad_fit.intercept;
0557 result.tbin_slope = tbin_fit.slope;
0558 result.tbin_intercept = tbin_fit.intercept;
0559
0560 if (fit_mode == Tpc_FittingTools::FIT_SAGITTA)
0561 {
0562 result.pad_sagitta = Tpc_FittingTools::fitSagitta(layer_pad_points);
0563 result.use_sagitta = result.pad_sagitta.ok;
0564 }
0565
0566 return result;
0567 }
0568
0569 FitResult fit_track_points(const std::vector<Tpc_ModuleTrackDisplay::HitPoint>& pts,
0570 const int fit_mode,
0571 const double weight_power,
0572 const double weight_floor_frac)
0573 {
0574 FitResult result;
0575 if (pts.size() < 2)
0576 {
0577 return result;
0578 }
0579
0580 double max_adc = 0.0;
0581 for (const auto& p : pts)
0582 {
0583 max_adc = std::max(max_adc, static_cast<double>(p.adc));
0584 }
0585 if (max_adc <= 0.0)
0586 {
0587 max_adc = 1.0;
0588 }
0589
0590 const double phi_ref = pts.front().phi;
0591 std::vector<Tpc_FittingTools::FitPoint> radius_phi_points;
0592 std::vector<Tpc_FittingTools::FitPoint> radius_tbin_points;
0593 radius_phi_points.reserve(pts.size());
0594 radius_tbin_points.reserve(pts.size());
0595
0596 result.rmin = pts.front().radius;
0597 result.rmax = pts.front().radius;
0598 result.mean_phi = 0.0;
0599
0600 for (const auto& p : pts)
0601 {
0602 if (!is_good_number(p.radius) || !is_good_number(p.phi))
0603 {
0604 continue;
0605 }
0606 const double w = Tpc_FittingTools::adcWeight(static_cast<double>(p.adc), max_adc,
0607 weight_power, weight_floor_frac);
0608 const double phi_unwrapped = unwrap_phi_near(p.phi, phi_ref);
0609 radius_phi_points.emplace_back(p.radius, phi_unwrapped, w);
0610 radius_tbin_points.emplace_back(p.radius, static_cast<double>(p.tbin), w);
0611 result.mean_phi += phi_unwrapped;
0612 result.rmin = std::min(result.rmin, p.radius);
0613 result.rmax = std::max(result.rmax, p.radius);
0614 }
0615
0616 if (radius_phi_points.empty())
0617 {
0618 return result;
0619 }
0620 result.mean_phi /= static_cast<double>(radius_phi_points.size());
0621
0622 const Tpc_FittingTools::LineFit phi_fit = Tpc_FittingTools::fitLine(radius_phi_points);
0623 const Tpc_FittingTools::LineFit tbin_fit = Tpc_FittingTools::fitLine(radius_tbin_points);
0624
0625 if (!phi_fit.ok || !tbin_fit.ok)
0626 {
0627 return result;
0628 }
0629
0630 result.ok = true;
0631 result.phi_slope = phi_fit.slope;
0632 result.phi_intercept = phi_fit.intercept;
0633 result.tbin_slope = tbin_fit.slope;
0634 result.tbin_intercept = tbin_fit.intercept;
0635
0636 if (fit_mode == Tpc_FittingTools::FIT_SAGITTA)
0637 {
0638 result.phi_sagitta = Tpc_FittingTools::fitSagitta(radius_phi_points);
0639 result.use_sagitta = result.phi_sagitta.ok;
0640 }
0641
0642 return result;
0643 }
0644
0645 void add_fit_objects(GraphBundle& b,
0646 const GroupKey& key,
0647 const unsigned int evt,
0648 const unsigned int tid,
0649 const int color,
0650 const FitResult& fit)
0651 {
0652 if (!fit.ok || fit.rmax <= fit.rmin)
0653 {
0654 return;
0655 }
0656
0657 const std::string tag = std::format("s{}_sec{:02}_mod{}_trk{}", key.side, key.sector, key.module, tid);
0658
0659 TGraph* g_phi_radius_fit = new TGraph();
0660 g_phi_radius_fit->SetName(std::format("g_{}_phi_radius_fit", tag).c_str());
0661 g_phi_radius_fit->SetTitle(std::format("event {} track {} fit;#phi;radius [cm]", evt, tid).c_str());
0662 style_fit_graph(g_phi_radius_fit, color);
0663
0664 TGraph* g_tbin_radius_fit = new TGraph();
0665 g_tbin_radius_fit->SetName(std::format("g_{}_tbin_radius_fit", tag).c_str());
0666 g_tbin_radius_fit->SetTitle(std::format("event {} track {} fit;timebin;radius [cm]", evt, tid).c_str());
0667 style_fit_graph(g_tbin_radius_fit, color);
0668
0669 TGraph* g_tbin_phi_fit = new TGraph();
0670 g_tbin_phi_fit->SetName(std::format("g_{}_tbin_phi_fit", tag).c_str());
0671 g_tbin_phi_fit->SetTitle(std::format("event {} track {} fit;timebin;#phi", evt, tid).c_str());
0672 style_fit_graph(g_tbin_phi_fit, color);
0673
0674 const int npts = 51;
0675 TPolyLine3D* line3 = new TPolyLine3D(npts);
0676 const std::string line_name = std::format("line3_{}_tbin_phi_radius_fit", tag);
0677
0678 for (int i = 0; i < npts; ++i)
0679 {
0680 const double f = static_cast<double>(i) / static_cast<double>(npts - 1);
0681 const double r = fit.rmin + f * (fit.rmax - fit.rmin);
0682 double phi = fit.phi_slope * r + fit.phi_intercept;
0683 if (fit.use_sagitta)
0684 {
0685 double phi_sagitta = 0.0;
0686 if (sagitta_y_at_x(r, fit.mean_phi, fit.phi_sagitta, phi_sagitta))
0687 {
0688 phi = phi_sagitta;
0689 }
0690 }
0691 const double tbin = fit.tbin_slope * r + fit.tbin_intercept;
0692 const double phi_draw = unwrap_phi_near(phi, b.phi_reference);
0693
0694 g_phi_radius_fit->SetPoint(i, phi_draw, r);
0695 g_tbin_radius_fit->SetPoint(i, tbin, r);
0696 g_tbin_phi_fit->SetPoint(i, tbin, phi_draw);
0697 line3->SetPoint(i, tbin, phi_draw, r);
0698 }
0699
0700 style_fit_line_3d(line3, color);
0701 b.fit_lines_3d.push_back(line3);
0702 b.fit_line_names_3d.push_back(line_name);
0703
0704 b.mg_phi_radius_fits->Add(g_phi_radius_fit, "L");
0705 b.mg_tbin_radius_fits->Add(g_tbin_radius_fit, "L");
0706 b.mg_tbin_phi_fits->Add(g_tbin_phi_fit, "L");
0707 }
0708
0709 void add_hardware_fit_objects(GraphBundle& b,
0710 const GroupKey& key,
0711 const unsigned int tid,
0712 const int color,
0713 const HardwareFitResult& fit)
0714 {
0715 if (!fit.ok || fit.layer_max <= fit.layer_min)
0716 {
0717 return;
0718 }
0719
0720 const std::string tag = std::format("s{}_sec{:02}_mod{}_trk{}", key.side, key.sector, key.module, tid);
0721
0722 const int npts = 51;
0723 TPolyLine3D* line3 = new TPolyLine3D(npts);
0724 const std::string line_name = std::format("line3_{}_layer_pad_tbin_fit", tag);
0725
0726 for (int i = 0; i < npts; ++i)
0727 {
0728 const double f = static_cast<double>(i) / static_cast<double>(npts - 1);
0729 const double layer = fit.layer_min + f * (fit.layer_max - fit.layer_min);
0730 double pad = fit.pad_slope * layer + fit.pad_intercept;
0731 if (fit.use_sagitta)
0732 {
0733 double pad_sagitta = 0.0;
0734 if (sagitta_y_at_x(layer, fit.mean_pad, fit.pad_sagitta, pad_sagitta))
0735 {
0736 pad = pad_sagitta;
0737 }
0738 }
0739 const double tbin = fit.tbin_slope * layer + fit.tbin_intercept;
0740 line3->SetPoint(i, tbin, pad, layer);
0741 }
0742
0743 style_fit_line_3d(line3, color);
0744 b.hardware_fit_lines_3d.push_back(line3);
0745 b.hardware_fit_line_names_3d.push_back(line_name);
0746 }
0747
0748 void write_bundle(GraphBundle& b, const GroupKey& key)
0749 {
0750 const std::string tag = std::format("s{}_sec{:02}_mod{}", key.side, key.sector, key.module);
0751
0752 if (b.h3_hardware_hits)
0753 {
0754 b.c3_hardware_hits_fits = new TCanvas(std::format("c3_{}_hardware_hits_fits", tag).c_str(),
0755 std::format("side {} sector {} module {} hardware hits and display fits",
0756 key.side, key.sector, key.module)
0757 .c_str(),
0758 1200, 900);
0759 b.h3_hardware_hits->Draw("BOX2Z");
0760 for (auto& i : b.hardware_fit_lines_3d)
0761 {
0762 if (!i)
0763 {
0764 continue;
0765 }
0766 i->Draw("same");
0767 }
0768 b.c3_hardware_hits_fits->Modified();
0769 b.c3_hardware_hits_fits->Update();
0770 b.c3_hardware_hits_fits->Write();
0771 }
0772
0773 if (b.h3_adc_hits)
0774 {
0775 b.c3_adc_hits_fits = new TCanvas(std::format("c3_{}_adc_hits_fits", tag).c_str(),
0776 std::format("side {} sector {} module {} ADC hits and display fits",
0777 key.side, key.sector, key.module)
0778 .c_str(),
0779 1200, 900);
0780 b.h3_adc_hits->Draw("BOX2Z");
0781 for (auto& i : b.fit_lines_3d)
0782 {
0783 if (!i)
0784 {
0785 continue;
0786 }
0787 i->Draw("same");
0788
0789 }
0790 b.c3_adc_hits_fits->Modified();
0791 b.c3_adc_hits_fits->Update();
0792 b.c3_adc_hits_fits->Write();
0793 }
0794
0795 if (b.h3_adc_unassociated_hits)
0796 {
0797 b.c3_adc_unassociated_hits = new TCanvas(std::format("c3_{}_adc_unassociated_hits", tag).c_str(),
0798 std::format("side {} sector {} module {} unassociated ADC hits",
0799 key.side, key.sector, key.module)
0800 .c_str(),
0801 1200, 900);
0802 b.h3_adc_unassociated_hits->Draw("BOX2Z");
0803 b.c3_adc_unassociated_hits->Modified();
0804 b.c3_adc_unassociated_hits->Update();
0805 b.c3_adc_unassociated_hits->Write();
0806 }
0807
0808
0809
0810
0811
0812
0813
0814
0815
0816
0817 }
0818 }
0819
0820 Tpc_ModuleTrackDisplay::HitPoint::HitPoint()
0821 : ok(false)
0822 , hitsetkey(0)
0823 , hitkey(0)
0824 , side(-1)
0825 , sector(0)
0826 , layer(0)
0827 , pad(0)
0828 , tbin(0)
0829 , adc(0)
0830 , radius(0.0)
0831 , phi(0.0)
0832 {
0833 }
0834
0835 Tpc_ModuleTrackDisplay::Tpc_ModuleTrackDisplay(const std::string& name,
0836 const std::string& outfilename,
0837 const std::string& trackNodeName,
0838 unsigned int maxEventDisplays)
0839 : SubsysReco(name)
0840 , m_outfilename(outfilename)
0841 , m_trackNodeName(trackNodeName)
0842 , m_maxEventDisplays(maxEventDisplays)
0843 , m_evt(0)
0844 , m_eventsSaved(0)
0845 , m_outfile(nullptr)
0846 , m_tracks(nullptr)
0847 , m_hits(nullptr)
0848 , m_idealPadMap(new IdealPadMap())
0849 , m_fitMode(Tpc_FittingTools::FIT_SAGITTA)
0850 , m_fitWeightPower(1.0)
0851 , m_fitWeightFloorFrac(0.05)
0852 {
0853 }
0854
0855 Tpc_ModuleTrackDisplay::~Tpc_ModuleTrackDisplay()
0856 {
0857 delete m_idealPadMap;
0858 m_idealPadMap = nullptr;
0859 }
0860
0861 int Tpc_ModuleTrackDisplay::Init(PHCompositeNode* )
0862 {
0863
0864 m_outfile = new TFile(m_outfilename.c_str(), "RECREATE");
0865 if (!m_outfile || m_outfile->IsZombie())
0866 {
0867 std::cerr << "Tpc_ModuleTrackDisplay::Init - cannot open output file "
0868 << m_outfilename << std::endl;
0869 return Fun4AllReturnCodes::ABORTRUN;
0870 }
0871
0872 if (!m_idealPadMap)
0873 {
0874 m_idealPadMap = new IdealPadMap();
0875 }
0876
0877 if (!m_idealPadMap->is_loaded() && m_idealPadMap->load_from_cdb(Verbosity()) != 0)
0878 {
0879 std::cerr << "Tpc_ModuleTrackDisplay::Init - failed to load IdealPadMap from CDB"
0880 << std::endl;
0881 return Fun4AllReturnCodes::ABORTRUN;
0882 }
0883
0884 m_outfile->mkdir("events");
0885 std::cout << "Tpc_ModuleTrackDisplay::Init - writing display fits to "
0886 << m_outfilename << std::endl;
0887
0888 return Fun4AllReturnCodes::EVENT_OK;
0889 }
0890
0891 int Tpc_ModuleTrackDisplay::process_event(PHCompositeNode* topNode)
0892 {
0893 ++m_evt;
0894
0895 if (!get_nodes(topNode))
0896 {
0897 return Fun4AllReturnCodes::EVENT_OK;
0898 }
0899
0900 const unsigned int ntracks = m_tracks ? m_tracks->size() : 0;
0901 if (m_eventsSaved >= m_maxEventDisplays)
0902 {
0903 return Fun4AllReturnCodes::EVENT_OK;
0904 }
0905
0906 TDirectory* eventsTop = m_outfile->GetDirectory("events");
0907 if (!eventsTop)
0908 {
0909 eventsTop = m_outfile->mkdir("events");
0910 }
0911
0912 eventsTop->cd();
0913 TDirectory* eventDir = eventsTop->mkdir(std::format("event_{:06}", m_evt).c_str());
0914 if (!eventDir)
0915 {
0916 std::cerr << "Tpc_ModuleTrackDisplay::process_event - failed to create event directory"
0917 << std::endl;
0918 return Fun4AllReturnCodes::EVENT_OK;
0919 }
0920 eventDir->cd();
0921
0922 std::map<GroupKey, GraphBundle> bundles;
0923 std::set<unsigned long long> associated_hit_ids;
0924
0925 for (unsigned int itrk = 0; itrk < ntracks; ++itrk)
0926 {
0927 const Tpc_ModuleTrack* trk = m_tracks->get_track(itrk);
0928 if (!trk)
0929 {
0930 continue;
0931 }
0932
0933 const unsigned int tid = trk->get_track_id();
0934 const int color = track_color(itrk);
0935
0936 std::map<GroupKey, std::vector<HitPoint> > points_by_group;
0937
0938 for (unsigned int ih = 0; ih < trk->size_hit_indices(); ++ih)
0939 {
0940 const Tpc_ModuleTrack::HitIndex idx = trk->get_hit_index(ih);
0941 const HitPoint p = make_hit_point(idx.first, idx.second);
0942 if (!p.ok)
0943 {
0944 continue;
0945 }
0946
0947 associated_hit_ids.insert(make_unique_hit_id(p.hitsetkey, p.hitkey));
0948
0949 const int module = layer_to_module(p.layer);
0950 if (module < 0)
0951 {
0952 continue;
0953 }
0954
0955 const GroupKey key(p.side, static_cast<int>(p.sector), module);
0956 points_by_group[key].push_back(p);
0957 }
0958
0959 for (auto& pg : points_by_group)
0960 {
0961 const GroupKey& key = pg.first;
0962 const std::vector<HitPoint>& pts = pg.second;
0963 if (pts.empty())
0964 {
0965 continue;
0966 }
0967
0968 GraphBundle& b = get_bundle(bundles, key, m_evt, m_idealPadMap);
0969
0970 TGraph* g_phi_radius_hits = new TGraph();
0971 g_phi_radius_hits->SetName(std::format("g_s{}_sec{:02}_mod{}_trk{}_phi_radius_hits",
0972 key.side, key.sector, key.module, tid)
0973 .c_str());
0974 g_phi_radius_hits->SetTitle(std::format("event {} track {} hits;#phi;radius [cm]", m_evt, tid).c_str());
0975 style_hit_graph(g_phi_radius_hits, color);
0976
0977 TGraph* g_tbin_radius_hits = new TGraph();
0978 g_tbin_radius_hits->SetName(std::format("g_s{}_sec{:02}_mod{}_trk{}_tbin_radius_hits",
0979 key.side, key.sector, key.module, tid)
0980 .c_str());
0981 g_tbin_radius_hits->SetTitle(std::format("event {} track {} hits;timebin;radius [cm]", m_evt, tid).c_str());
0982 style_hit_graph(g_tbin_radius_hits, color);
0983
0984 TGraph* g_tbin_phi_hits = new TGraph();
0985 g_tbin_phi_hits->SetName(std::format("g_s{}_sec{:02}_mod{}_trk{}_tbin_phi_hits",
0986 key.side, key.sector, key.module, tid)
0987 .c_str());
0988 g_tbin_phi_hits->SetTitle(std::format("event {} track {} hits;timebin;#phi", m_evt, tid).c_str());
0989 style_hit_graph(g_tbin_phi_hits, color);
0990
0991 for (const auto& p : pts)
0992 {
0993 const int n = g_phi_radius_hits->GetN();
0994 const double phi_draw = unwrap_phi_near(p.phi, b.phi_reference);
0995
0996 g_phi_radius_hits->SetPoint(n, phi_draw, p.radius);
0997 g_tbin_radius_hits->SetPoint(n, static_cast<double>(p.tbin), p.radius);
0998 g_tbin_phi_hits->SetPoint(n, static_cast<double>(p.tbin), phi_draw);
0999
1000 const unsigned long long uid = make_unique_hit_id(p.hitsetkey, p.hitkey);
1001 if (b.filled_hit_ids.insert(uid).second)
1002 {
1003 if (b.h3_hardware_hits)
1004 {
1005 b.h3_hardware_hits->Fill(static_cast<double>(p.tbin),
1006 static_cast<double>(p.pad),
1007 static_cast<double>(p.layer),
1008 static_cast<double>(p.adc));
1009 }
1010 if (b.h3_adc_hits)
1011 {
1012 b.h3_adc_hits->Fill(static_cast<double>(p.tbin),
1013 phi_draw,
1014 p.radius,
1015 static_cast<double>(p.adc));
1016 }
1017 }
1018 }
1019
1020 b.mg_phi_radius_hits->Add(g_phi_radius_hits, "LP");
1021 b.mg_tbin_radius_hits->Add(g_tbin_radius_hits, "LP");
1022 b.mg_tbin_phi_hits->Add(g_tbin_phi_hits, "LP");
1023
1024 const FitResult fit = fit_track_points(pts, m_fitMode, m_fitWeightPower, m_fitWeightFloorFrac);
1025 add_fit_objects(b, key, m_evt, tid, color, fit);
1026
1027 const HardwareFitResult hardwareFit =
1028 fit_hardware_track_points(pts, m_fitMode, m_fitWeightPower, m_fitWeightFloorFrac);
1029 add_hardware_fit_objects(b, key, tid, color, hardwareFit);
1030 }
1031 }
1032
1033
1034 TrkrHitSetContainer::ConstRange hitset_range = m_hits->getHitSets();
1035 for (TrkrHitSetContainer::ConstIterator hsiter = hitset_range.first;
1036 hsiter != hitset_range.second; ++hsiter)
1037 {
1038 TrkrHitSet* hitset = hsiter->second;
1039 if (!hitset)
1040 {
1041 continue;
1042 }
1043
1044 const TrkrDefs::hitsetkey hsk = hsiter->first;
1045 const unsigned int layer = TrkrDefs::getLayer(hsk);
1046 const int module = layer_to_module(layer);
1047 if (module < 0)
1048 {
1049 continue;
1050 }
1051
1052 const int side = static_cast<int>(TpcDefs::getSide(hsk));
1053 const int sector = static_cast<int>(TpcDefs::getSectorId(hsk));
1054 const GroupKey key(side, sector, module);
1055 GraphBundle& b = get_bundle(bundles, key, m_evt, m_idealPadMap);
1056
1057 TrkrHitSet::ConstRange hit_range = hitset->getHits();
1058 for (TrkrHitSet::ConstIterator hiter = hit_range.first;
1059 hiter != hit_range.second; ++hiter)
1060 {
1061 const TrkrDefs::hitkey hk = hiter->first;
1062 const unsigned long long uid = make_unique_hit_id(hsk, hk);
1063 if (associated_hit_ids.contains(uid))
1064 {
1065 continue;
1066 }
1067
1068 const HitPoint p = make_hit_point(hsk, hk);
1069 if (!p.ok)
1070 {
1071 continue;
1072 }
1073
1074 if (b.h3_adc_unassociated_hits)
1075 {
1076 b.h3_adc_unassociated_hits->Fill(static_cast<double>(p.tbin),
1077 unwrap_phi_near(p.phi, b.phi_reference),
1078 p.radius,
1079 static_cast<double>(p.adc));
1080 }
1081 }
1082 }
1083
1084 for (auto& bundle : bundles)
1085 {
1086 eventDir->cd();
1087 write_bundle(bundle.second, bundle.first);
1088 }
1089
1090 std::cout << "Tpc_ModuleTrackDisplay - saved event "
1091 << m_evt << " with " << ntracks
1092 << " tracks and " << bundles.size()
1093 << " side/sector/module groups" << std::endl;
1094
1095 ++m_eventsSaved;
1096 return Fun4AllReturnCodes::EVENT_OK;
1097 }
1098
1099 int Tpc_ModuleTrackDisplay::End(PHCompositeNode* )
1100 {
1101 if (m_outfile)
1102 {
1103 std::cout << "Tpc_ModuleTrackDisplay::End - events seen: " << m_evt << std::endl;
1104 m_outfile->cd();
1105 m_outfile->Write();
1106 m_outfile->Close();
1107 delete m_outfile;
1108 m_outfile = nullptr;
1109 }
1110
1111 std::cout << "Tpc_ModuleTrackDisplay::End - wrote " << m_outfilename << std::endl;
1112 return Fun4AllReturnCodes::EVENT_OK;
1113 }
1114
1115 bool Tpc_ModuleTrackDisplay::get_nodes(PHCompositeNode* topNode)
1116 {
1117 m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1118 m_tracks = findNode::getClass<Tpc_ModuleTrackContainer>(topNode, m_trackNodeName);
1119
1120 if (!m_tracks)
1121 {
1122 const char* candidate_names[] = {
1123 "TPC_MODULETRACKS",
1124 "Tpc_ModuleTrackReco",
1125 "Tpc_ModuleTrackContainer",
1126 "TPC_MODULETRACKCONTAINER",
1127 "TPC_MODULETRACKS_CONTAINER"};
1128
1129 for (unsigned int i = 0;
1130 i < sizeof(candidate_names) / sizeof(candidate_names[0]) && !m_tracks;
1131 ++i)
1132 {
1133 m_tracks = findNode::getClass<Tpc_ModuleTrackContainer>(topNode, candidate_names[i]);
1134 if (m_tracks)
1135 {
1136 m_trackNodeName = candidate_names[i];
1137 }
1138 }
1139 }
1140
1141
1142
1143
1144 if (!m_idealPadMap)
1145 {
1146 m_idealPadMap = new IdealPadMap();
1147 }
1148
1149 if (!m_tracks)
1150 {
1151 std::cerr << "Tpc_ModuleTrackDisplay - could not find Tpc_ModuleTrackContainer node" << std::endl;
1152 return false;
1153 }
1154
1155 if (!m_hits)
1156 {
1157 std::cerr << "Tpc_ModuleTrackDisplay - missing TRKR_HITSET" << std::endl;
1158 return false;
1159 }
1160
1161 if (!m_idealPadMap)
1162 {
1163 std::cerr << "Tpc_ModuleTrackDisplay - failed to create IdealPadMap" << std::endl;
1164 return false;
1165 }
1166
1167 return true;
1168 }
1169
1170 Tpc_ModuleTrackDisplay::HitPoint
1171 Tpc_ModuleTrackDisplay::make_hit_point(const TrkrDefs::hitsetkey hsk,
1172 const TrkrDefs::hitkey hk) const
1173 {
1174 HitPoint p;
1175
1176 TrkrHitSet* hitset = m_hits ? m_hits->findHitSet(hsk) : nullptr;
1177 if (!hitset)
1178 {
1179 return p;
1180 }
1181
1182 TrkrHit* hit = hitset->getHit(hk);
1183 if (!hit)
1184 {
1185 return p;
1186 }
1187
1188 p.hitsetkey = hsk;
1189 p.hitkey = hk;
1190 p.side = static_cast<int>(TpcDefs::getSide(hsk));
1191 p.sector = static_cast<unsigned int>(TpcDefs::getSectorId(hsk));
1192 p.layer = TrkrDefs::getLayer(hsk);
1193 p.pad = TpcDefs::getPad(hk);
1194 p.tbin = TpcDefs::getTBin(hk);
1195 p.adc = hit->getAdc();
1196
1197 if (layer_to_module(p.layer) < 0)
1198 {
1199 return p;
1200 }
1201
1202 p.radius = ideal_radius(p.layer);
1203 p.phi = ideal_phi(p.side, p.sector, p.layer, p.pad);
1204
1205 if (!is_good_number(p.radius) || !is_good_number(p.phi))
1206 {
1207 return p;
1208 }
1209
1210 p.ok = true;
1211 return p;
1212 }
1213
1214 double Tpc_ModuleTrackDisplay::ideal_radius(const unsigned int layer) const
1215 {
1216 if (!m_idealPadMap)
1217 {
1218 return 0.0;
1219 }
1220
1221
1222
1223 return m_idealPadMap->get_radius(layer);
1224 }
1225
1226 double Tpc_ModuleTrackDisplay::ideal_phi(const int side,
1227 const unsigned int sector,
1228 const unsigned int layer,
1229 const unsigned int pad) const
1230 {
1231 if (!m_idealPadMap)
1232 {
1233 return 0.0;
1234 }
1235
1236
1237
1238 const unsigned int pads_per_sector = m_idealPadMap->get_pads_per_sector_for_layer(layer);
1239 if (pads_per_sector == 0U)
1240 {
1241 return 0.0;
1242 }
1243 const unsigned int local_pad = pad % pads_per_sector;
1244 return m_idealPadMap->get_phi(side, sector, layer, local_pad);
1245 }