Back to home page

sPhenix code displayed by LXR

 
 

    


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   // NOLINTBEGIN(misc-non-private-member-variables-in-classes)
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     // x = timebin, y = pad, z = layer, weight = ADC
0340     TH3D* h3_hardware_hits;
0341 
0342     // x = timebin, y = phi, z = radius, weight = ADC
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   // NOLINTEND(misc-non-private-member-variables-in-classes)
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     // first and last pad centers in this sector
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);  // NOLINT(readability-suspicious-call-argument)
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         // b.fit_lines_3d[i]->Write(b.fit_line_names_3d[i].c_str());
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         if (b.mg_phi_radius_hits) b.mg_phi_radius_hits->Write();
0809         if (b.mg_tbin_radius_hits) b.mg_tbin_radius_hits->Write();
0810         if (b.mg_tbin_phi_hits) b.mg_tbin_phi_hits->Write();
0811         if (b.mg_phi_radius_fits) b.mg_phi_radius_fits->Write();
0812         if (b.mg_tbin_radius_fits) b.mg_tbin_radius_fits->Write();
0813         if (b.mg_tbin_phi_fits) b.mg_tbin_phi_fits->Write();
0814         if (b.h3_adc_hits) b.h3_adc_hits->Write();
0815         if (b.h3_adc_unassociated_hits) b.h3_adc_unassociated_hits->Write();
0816         */
0817   }
0818 }  // namespace
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* /*unused*/)
0862 {
0863   // cppcheck-suppress publicAllocationError
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   // Fill all TRKR_HITSET hits that are not associated with any reconstructed in-module track.
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* /*unused*/)
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   // IdealPadMap is not a Fun4All node/container here.
1142   // The display owns it and uses it only to translate hardware coordinates
1143   // (side, sector, layer, pad) -> (phi, radius).
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   // IdealPadMap is intentionally the only geometry source used by this display.
1222   // If your IdealPadMap method names differ, this is the only radius call to adjust.
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   // Keep HitPoint::pad in raw hardware coordinates for display, but IdealPadMap
1237   // expects the pad index local to the sector.
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 }