Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-08-31 08:21:22 

0001 #include "Tpc_ModuleTrackReco.h"
0002 #include "Tpc_FittingTools.h"
0003 
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/PHIODataNode.h>
0008 #include <phool/PHNodeIterator.h>
0009 #include <phool/PHObject.h>
0010 #include <phool/getClass.h>
0011 
0012 #include <TFile.h>
0013 #include <TTree.h>
0014 
0015 #include "Tpc_ModuleTrack.h"
0016 #include "Tpc_ModuleTrackContainer.h"
0017 #include "Tpc_ModuleTrackContainerv1.h"
0018 #include "Tpc_ModuleTrackv1.h"
0019 
0020 #include <trackbase/TpcDefs.h>
0021 #include <trackbase/TrkrDefs.h>
0022 #include <trackbase/TrkrHit.h>
0023 #include <trackbase/TrkrHitSet.h>
0024 #include <trackbase/TrkrHitSetContainer.h>
0025 
0026 #include <pthread.h>
0027 
0028 #include <algorithm>
0029 #include <cmath>
0030 #include <cstdint>
0031 #include <deque>
0032 #include <iostream>
0033 #include <limits>
0034 #include <map>
0035 #include <utility>
0036 #include <vector>
0037 
0038 // ===================================================================
0039 // Internal helpers (anonymous namespace)
0040 // ===================================================================
0041 namespace
0042 {
0043   // NOLINTBEGIN(misc-non-private-member-variables-in-classes)
0044   struct BlobAdcSort
0045   {
0046     const std::vector<InModuleThreadData::Blob>* blobs;
0047     explicit BlobAdcSort(const std::vector<InModuleThreadData::Blob>* b)
0048       : blobs(b)
0049     {
0050     }
0051     bool operator()(unsigned int a, unsigned int b) const
0052     {
0053       return (*blobs)[a].adc > (*blobs)[b].adc;
0054     }
0055   };
0056 
0057   // Generic fitting functions live in Tpc_FittingTools.{h,cc}.
0058 
0059   bool fit_track_from_blobs(const std::vector<InModuleThreadData::Blob>& blobs,
0060                             const std::vector<unsigned int>& idx,
0061                             double weight_power,
0062                             double floor_frac,
0063                             InModuleThreadData::Track& trk)
0064   {
0065     if (idx.size() < 2)
0066     {
0067       return false;
0068     }
0069 
0070     double maxadc = 0.0;
0071     for (unsigned int i : idx)
0072     {
0073       maxadc = std::max(blobs[i].adc, maxadc);
0074     }
0075 
0076     std::vector<double> x;
0077     std::vector<double> pad;
0078     std::vector<double> tbin;
0079     std::vector<double> w;
0080     x.reserve(idx.size());
0081     pad.reserve(idx.size());
0082     tbin.reserve(idx.size());
0083     w.reserve(idx.size());
0084 
0085     unsigned int first_layer = 999999;
0086     unsigned int last_layer = 0;
0087 
0088     for (unsigned int i : idx)
0089     {
0090       const InModuleThreadData::Blob& bl = blobs[i];
0091       x.push_back(static_cast<double>(bl.layer));
0092       pad.push_back(bl.pad);
0093       tbin.push_back(bl.tbin);
0094       w.push_back(Tpc_FittingTools::adcWeight(bl.adc, maxadc, weight_power, floor_frac));
0095 
0096       first_layer = std::min(bl.layer, first_layer);
0097       last_layer = std::max(bl.layer, last_layer);
0098     }
0099 
0100     double mp = 0.0;
0101     double bp = 0.0;
0102     double cp = 0.0;
0103     double mt = 0.0;
0104     double bt = 0.0;
0105     double ct = 0.0;
0106     int ndp = 0;
0107     int ndt = 0;
0108 
0109     if (!Tpc_FittingTools::weightedLineFit(x, pad, w, mp, bp, cp, ndp))
0110     {
0111       return false;
0112     }
0113     if (!Tpc_FittingTools::weightedLineFit(x, tbin, w, mt, bt, ct, ndt))
0114     {
0115       return false;
0116     }
0117 
0118     trk.first_layer = first_layer;
0119     trk.last_layer = last_layer;
0120     trk.nblobs = static_cast<unsigned int>(idx.size());
0121     trk.nrawhits = 0;
0122     trk.pad_slope = mp;
0123     trk.pad_intercept = bp;
0124     trk.tbin_slope = mt;
0125     trk.tbin_intercept = bt;
0126     trk.blob_indices = idx;
0127     trk.raw_hit_indices.clear();
0128 
0129     return true;
0130   }
0131 
0132   void collect_raw_indices_from_blob_chain(const std::vector<InModuleThreadData::Blob>& blobs,
0133                                            const std::vector<unsigned int>& blob_idx,
0134                                            std::vector<unsigned int>& raw_idx)
0135   {
0136     raw_idx.clear();
0137     for (unsigned int ib : blob_idx)
0138     {
0139       const InModuleThreadData::Blob& bl = blobs[ib];
0140       for (unsigned int raw_hit_indice : bl.raw_hit_indices)
0141       {
0142         raw_idx.push_back(raw_hit_indice);
0143       }
0144     }
0145   }
0146 
0147   bool make_track_from_blob_chain(const std::vector<InModuleThreadData::RawHit>& raw_hits,
0148                                   const std::vector<InModuleThreadData::Blob>& blobs,
0149                                   const std::vector<unsigned int>& blob_idx,
0150                                   double weight_power,
0151                                   double floor_frac,
0152                                   InModuleThreadData::Track& trk)
0153   {
0154     if (!fit_track_from_blobs(blobs, blob_idx, weight_power, floor_frac, trk))
0155     {
0156       return false;
0157     }
0158 
0159     std::vector<unsigned int> raw_idx;
0160     collect_raw_indices_from_blob_chain(blobs, blob_idx, raw_idx);
0161     if (raw_idx.empty())
0162     {
0163       return false;
0164     }
0165 
0166     unsigned int first_layer = 999999;
0167     unsigned int last_layer = 0;
0168     for (unsigned int i : raw_idx)
0169     {
0170       const InModuleThreadData::RawHit& rh = raw_hits[i];
0171       first_layer = std::min(rh.layer, first_layer);
0172       last_layer = std::max(rh.layer, last_layer);
0173     }
0174 
0175     trk.first_layer = first_layer;
0176     trk.last_layer = last_layer;
0177     trk.nblobs = static_cast<unsigned int>(blob_idx.size());
0178     trk.nrawhits = static_cast<unsigned int>(raw_idx.size());
0179     trk.blob_indices = blob_idx;
0180     trk.raw_hit_indices = raw_idx;
0181 
0182     return true;
0183   }
0184 
0185   // -------------------------------------------------------------------
0186   // Track-piece connection
0187   // -------------------------------------------------------------------
0188   struct TrackStartSort
0189   {
0190     const std::vector<InModuleThreadData::Track>* tracks;
0191     explicit TrackStartSort(const std::vector<InModuleThreadData::Track>* t)
0192       : tracks(t)
0193     {
0194     }
0195     bool operator()(unsigned int a, unsigned int b) const
0196     {
0197       const InModuleThreadData::Track& ta = (*tracks)[a];
0198       const InModuleThreadData::Track& tb = (*tracks)[b];
0199       if (ta.first_layer != tb.first_layer)
0200       {
0201         return ta.first_layer < tb.first_layer;
0202       }
0203       if (ta.last_layer != tb.last_layer)
0204       {
0205         return ta.last_layer < tb.last_layer;
0206       }
0207       return ta.nrawhits > tb.nrawhits;
0208     }
0209   };
0210 
0211   void append_unique_blob_indices(std::vector<unsigned int>& dst,
0212                                   const std::vector<unsigned int>& src)
0213   {
0214     for (unsigned int i : src)
0215     {
0216       if (std::find(dst.begin(), dst.end(), i) == dst.end())
0217       {
0218         dst.push_back(i);
0219       }
0220     }
0221   }
0222 
0223   bool tracks_can_connect(const InModuleThreadData::Track& a,
0224                           const InModuleThreadData::Track& b,
0225                           unsigned int connect_max_layer_gap,
0226                           double connect_dp,
0227                           double connect_dt,
0228                           double connect_dpad_slope,
0229                           double connect_dtbin_slope,
0230                           double& score)
0231   {
0232     score = std::numeric_limits<double>::max();
0233 
0234     if (a.last_layer >= b.first_layer)
0235     {
0236       return false;
0237     }
0238 
0239     const unsigned int gap = b.first_layer - a.last_layer - 1;
0240     if (gap > connect_max_layer_gap)
0241     {
0242       return false;
0243     }
0244 
0245     const double lmatch = 0.5 * (static_cast<double>(a.last_layer) +
0246                                  static_cast<double>(b.first_layer));
0247 
0248     const double pad_a = a.pad_slope * lmatch + a.pad_intercept;
0249     const double pad_b = b.pad_slope * lmatch + b.pad_intercept;
0250     const double tbin_a = a.tbin_slope * lmatch + a.tbin_intercept;
0251     const double tbin_b = b.tbin_slope * lmatch + b.tbin_intercept;
0252 
0253     const double dp = std::fabs(pad_a - pad_b);
0254     const double dt = std::fabs(tbin_a - tbin_b);
0255     const double dmp = std::fabs(a.pad_slope - b.pad_slope);
0256     const double dmt = std::fabs(a.tbin_slope - b.tbin_slope);
0257 
0258     if (dp > connect_dp)
0259     {
0260       return false;
0261     }
0262     if (dt > connect_dt)
0263     {
0264       return false;
0265     }
0266     if (dmp > connect_dpad_slope)
0267     {
0268       return false;
0269     }
0270     if (dmt > connect_dtbin_slope)
0271     {
0272       return false;
0273     }
0274 
0275     score = (dp / connect_dp) * (dp / connect_dp) + (dt / connect_dt) * (dt / connect_dt) + (dmp / connect_dpad_slope) * (dmp / connect_dpad_slope) + (dmt / connect_dtbin_slope) * (dmt / connect_dtbin_slope) + 0.05 * static_cast<double>(gap);
0276 
0277     return true;
0278   }
0279 
0280   void connect_track_pieces_in_module(InModuleThreadData* d)
0281   {
0282     if (!d || d->tracks.size() < 2)
0283     {
0284       return;
0285     }
0286 
0287     std::vector<InModuleThreadData::Track> pieces = d->tracks;
0288     std::vector<InModuleThreadData::Track> output;
0289     std::vector<int> used(pieces.size(), 0);
0290 
0291     std::vector<unsigned int> order;
0292     order.reserve(pieces.size());
0293     for (unsigned int i = 0; i < pieces.size(); ++i)
0294     {
0295       order.push_back(i);
0296     }
0297     std::sort(order.begin(), order.end(), TrackStartSort(&pieces));
0298 
0299     for (unsigned int io = 0; io < order.size(); ++io)
0300     {
0301       const unsigned int iseed = order[io];
0302       if (used[iseed])
0303       {
0304         continue;
0305       }
0306 
0307       InModuleThreadData::Track current = pieces[iseed];
0308       used[iseed] = 1;
0309 
0310       bool merged_any = true;
0311       while (merged_any)
0312       {
0313         merged_any = false;
0314         int best_j = -1;
0315         double best_score = std::numeric_limits<double>::max();
0316 
0317         for (unsigned int j : order)
0318         {
0319           if (used[j])
0320           {
0321             continue;
0322           }
0323 
0324           double score = 0.0;
0325           if (!tracks_can_connect(current, pieces[j],
0326                                   d->connect_max_layer_gap,
0327                                   d->connect_dp,
0328                                   d->connect_dt,
0329                                   d->connect_dpad_slope,
0330                                   d->connect_dtbin_slope,
0331                                   score))
0332           {
0333             continue;
0334           }
0335 
0336           if (score < best_score)
0337           {
0338             best_score = score;
0339             best_j = static_cast<int>(j);
0340           }
0341         }
0342 
0343         if (best_j >= 0)
0344         {
0345           append_unique_blob_indices(current.blob_indices, pieces[best_j].blob_indices);
0346 
0347           InModuleThreadData::Track refit;
0348           if (make_track_from_blob_chain(d->raw_hits, d->blobs,
0349                                          current.blob_indices,
0350                                          d->weight_power,
0351                                          d->adc_weight_floor_frac,
0352                                          refit))
0353           {
0354             current = refit;
0355             used[best_j] = 1;
0356             merged_any = true;
0357           }
0358         }
0359       }
0360 
0361       output.push_back(current);
0362     }
0363 
0364     for (unsigned int i = 0; i < output.size(); ++i)
0365     {
0366       output[i].track_id = i;
0367     }
0368 
0369     if (d->verbosity > 1)
0370     {
0371       std::cout << "Tpc_ModuleTrackReco connect pieces: region=" << d->region
0372                 << " sector=" << d->sector << " side=" << d->side
0373                 << " pieces=" << pieces.size()
0374                 << " connected_tracks=" << output.size() << std::endl;
0375     }
0376 
0377     d->tracks.swap(output);
0378   }
0379 
0380   // -------------------------------------------------------------------
0381   // Per-module worker functions (exact originals)
0382   // -------------------------------------------------------------------
0383 
0384   struct RawHitTimeSort
0385   {
0386     const std::vector<InModuleThreadData::RawHit>* raw_hits;
0387     explicit RawHitTimeSort(const std::vector<InModuleThreadData::RawHit>* h)
0388       : raw_hits(h)
0389     {
0390     }
0391     bool operator()(unsigned int a, unsigned int b) const
0392     {
0393       const InModuleThreadData::RawHit& ha = (*raw_hits)[a];
0394       const InModuleThreadData::RawHit& hb = (*raw_hits)[b];
0395       if (ha.tbin != hb.tbin)
0396       {
0397         return ha.tbin < hb.tbin;
0398       }
0399       return ha.adc > hb.adc;
0400     }
0401   };
0402   // NOLINTEND(misc-non-private-member-variables-in-classes)
0403 
0404   void reject_long_pad_noise(InModuleThreadData* d)
0405   {
0406     if (!d)
0407     {
0408       return;
0409     }
0410     if (d->noise_max_consecutive_timebins <= 0)
0411     {
0412       return;
0413     }
0414     d->noise_keep_first_timebins = std::max(d->noise_keep_first_timebins, 0);
0415     if (d->raw_hits.empty())
0416     {
0417       return;
0418     }
0419 
0420     using LayerPadKey = std::pair<unsigned int, unsigned short>;
0421     std::map<LayerPadKey, std::vector<unsigned int> > by_layer_pad;
0422 
0423     for (unsigned int i = 0; i < d->raw_hits.size(); ++i)
0424     {
0425       const InModuleThreadData::RawHit& rh = d->raw_hits[i];
0426       by_layer_pad[LayerPadKey(rh.layer, rh.pad)].push_back(i);
0427     }
0428 
0429     std::vector<int> remove(d->raw_hits.size(), 0);
0430     unsigned int nremoved = 0;
0431 
0432     for (auto& it : by_layer_pad)
0433     {
0434       std::vector<unsigned int>& idx = it.second;
0435       if (idx.size() <= static_cast<unsigned int>(d->noise_max_consecutive_timebins))
0436       {
0437         continue;
0438       }
0439 
0440       std::sort(idx.begin(), idx.end(), RawHitTimeSort(&d->raw_hits));
0441 
0442       unsigned int run_start = 0;
0443       while (run_start < idx.size())
0444       {
0445         unsigned int run_end = run_start;
0446         while (run_end + 1 < idx.size())
0447         {
0448           const unsigned short t0 = d->raw_hits[idx[run_end]].tbin;
0449           const unsigned short t1 = d->raw_hits[idx[run_end + 1]].tbin;
0450           if (static_cast<int>(t1) != static_cast<int>(t0) + 1)
0451           {
0452             break;
0453           }
0454           ++run_end;
0455         }
0456 
0457         const unsigned int run_len = run_end - run_start + 1;
0458         if (run_len > static_cast<unsigned int>(d->noise_max_consecutive_timebins))
0459         {
0460           const unsigned int keep_until =
0461               run_start + static_cast<unsigned int>(d->noise_keep_first_timebins);
0462 
0463           // bool in_non_increasing_tail = true;
0464           for (unsigned int ir = run_start; ir <= run_end; ++ir)
0465           {
0466             if (ir < keep_until)
0467             {
0468               continue;
0469             }
0470 
0471             const unsigned int cur_idx = idx[ir];
0472             const unsigned int prev_idx = idx[ir - 1];
0473 
0474             const unsigned short cur_adc = d->raw_hits[cur_idx].adc;
0475             const unsigned short prev_adc = d->raw_hits[prev_idx].adc;
0476 
0477             const int cur_adc_i = static_cast<int>(cur_adc);
0478             const int prev_adc_i = static_cast<int>(prev_adc);
0479             const int tol = d->noise_adc_tolerance;
0480 
0481             // remove tail if ADC is flat, decreasing, or only mildly increasing
0482             if (cur_adc_i <= prev_adc_i + tol)
0483             {
0484               if (!remove[cur_idx])
0485               {
0486                 remove[cur_idx] = 1;
0487                 ++nremoved;
0488               }
0489             }
0490             else
0491             {
0492               // Preserve a later rising signal on top of the long pad tail.
0493               // in_non_increasing_tail = false;
0494             }
0495           }
0496         }
0497 
0498         run_start = run_end + 1;
0499       }
0500     }
0501 
0502     if (nremoved == 0)
0503     {
0504       return;
0505     }
0506 
0507     std::vector<InModuleThreadData::RawHit> kept;
0508     kept.reserve(d->raw_hits.size() - nremoved);
0509     for (unsigned int i = 0; i < d->raw_hits.size(); ++i)
0510     {
0511       if (!remove[i])
0512       {
0513         kept.push_back(d->raw_hits[i]);
0514       }
0515     }
0516 
0517     d->raw_hits.swap(kept);
0518 
0519     if (d->verbosity > 1)
0520     {
0521       std::cout << "Tpc_ModuleTrackReco noise rejection: region=" << d->region
0522                 << " sector=" << d->sector << " side=" << d->side
0523                 << " removed " << nremoved << " long same-pad tail hits"
0524                 << std::endl;
0525     }
0526   }
0527 
0528   void collect_raw_hits(InModuleThreadData* d)
0529   {
0530     d->raw_hits.clear();
0531 
0532     for (unsigned int ihs = 0; ihs < d->layer_hitsets.size(); ++ihs)
0533     {
0534       TrkrHitSet* hitset = d->layer_hitsets[ihs].hitset;
0535       if (!hitset)
0536       {
0537         continue;
0538       }
0539 
0540       TrkrHitSet::ConstRange range = hitset->getHits();
0541       for (TrkrHitSet::ConstIterator hitr = range.first; hitr != range.second; ++hitr)
0542       {
0543         const TrkrDefs::hitkey hitkey = hitr->first;
0544         const TrkrHit* hit = hitr->second;
0545         if (!hit)
0546         {
0547           continue;
0548         }
0549 
0550         const unsigned short pad = TpcDefs::getPad(hitkey);
0551         const unsigned short tbin = TpcDefs::getTBin(hitkey);
0552         const unsigned short rawAdc = hit->getAdc();
0553         const double fadc = static_cast<double>(rawAdc) - d->pedestal;
0554         if (fadc <= 0.0)
0555         {
0556           continue;
0557         }
0558 
0559         InModuleThreadData::RawHit rh;
0560         rh.layer = d->layer_hitsets[ihs].layer;
0561         rh.hitsetkey = d->layer_hitsets[ihs].hitsetkey;
0562         rh.hitkey = hitkey;
0563         rh.pad = pad;
0564         rh.tbin = tbin;
0565         rh.adc = static_cast<unsigned short>(fadc);
0566         d->raw_hits.push_back(rh);
0567       }
0568     }
0569   }
0570 
0571   void build_blobs(InModuleThreadData* d)
0572   {
0573     d->blobs.clear();
0574     const unsigned int n = static_cast<unsigned int>(d->raw_hits.size());
0575     std::vector<int> used(n, 0);
0576 
0577     for (unsigned int i = 0; i < n; ++i)
0578     {
0579       if (used[i])
0580       {
0581         continue;
0582       }
0583 
0584       used[i] = 1;
0585       std::deque<unsigned int> q;
0586       q.push_back(i);
0587 
0588       double sw = 0.0;
0589       double sp = 0.0;
0590       double st = 0.0;
0591       unsigned int nh = 0;
0592       const unsigned int layer = d->raw_hits[i].layer;
0593 
0594       InModuleThreadData::Blob bl;
0595 
0596       while (!q.empty())
0597       {
0598         const unsigned int a = q.front();
0599         q.pop_front();
0600 
0601         const InModuleThreadData::RawHit& ha = d->raw_hits[a];
0602         bl.raw_hit_indices.push_back(a);
0603         const double wa = static_cast<double>(ha.adc);
0604         sw += wa;
0605         sp += wa * static_cast<double>(ha.pad);
0606         st += wa * static_cast<double>(ha.tbin);
0607         ++nh;
0608 
0609         for (unsigned int j = 0; j < n; ++j)
0610         {
0611           if (used[j])
0612           {
0613             continue;
0614           }
0615           const InModuleThreadData::RawHit& hb = d->raw_hits[j];
0616           if (hb.layer != layer)
0617           {
0618             continue;
0619           }
0620 
0621           const int dp = std::abs(static_cast<int>(hb.pad) - static_cast<int>(ha.pad));
0622           const int dt = std::abs(static_cast<int>(hb.tbin) - static_cast<int>(ha.tbin));
0623           if (dp <= d->blob_dp && dt <= d->blob_dt)
0624           {
0625             used[j] = 1;
0626             q.push_back(j);
0627           }
0628         }
0629       }
0630 
0631       if (sw <= 0.0)
0632       {
0633         continue;
0634       }
0635 
0636       bl.layer = layer;
0637       bl.pad = sp / sw;
0638       bl.tbin = st / sw;
0639       bl.adc = sw;
0640       bl.nhits = nh;
0641       bl.used = 0;
0642       d->blobs.push_back(bl);
0643     }
0644   }
0645 
0646   int find_best_blob_on_layer(const InModuleThreadData* d,
0647                               unsigned int target_layer,
0648                               double pred_pad,
0649                               double pred_tbin)
0650   {
0651     int best = -1;
0652     double best_score = std::numeric_limits<double>::max();
0653 
0654     for (unsigned int i = 0; i < d->blobs.size(); ++i)
0655     {
0656       const InModuleThreadData::Blob& bl = d->blobs[i];
0657       if (bl.used)
0658       {
0659         continue;
0660       }
0661       if (bl.layer != target_layer)
0662       {
0663         continue;
0664       }
0665 
0666       const double dp = bl.pad - pred_pad;
0667       const double dt = bl.tbin - pred_tbin;
0668       if (std::fabs(dp) > d->search_dp)
0669       {
0670         continue;
0671       }
0672       if (std::fabs(dt) > d->search_dt)
0673       {
0674         continue;
0675       }
0676 
0677       const double score = (dp * dp) / (d->search_dp * d->search_dp + 1.0e-9) + (dt * dt) / (d->search_dt * d->search_dt + 1.0e-9) - 0.01 * std::log(bl.adc + 1.0);
0678       if (score < best_score)
0679       {
0680         best_score = score;
0681         best = static_cast<int>(i);
0682       }
0683     }
0684 
0685     return best;
0686   }
0687 
0688   void grow_one_direction(InModuleThreadData* d,
0689                           std::vector<unsigned int>& chain,
0690                           int direction)
0691   {
0692     while (true)
0693     {
0694       unsigned int edge_layer = d->blobs[chain.back()].layer;
0695       if (direction < 0)
0696       {
0697         edge_layer = d->blobs[chain.front()].layer;
0698       }
0699 
0700       if (direction > 0 && edge_layer >= 54)
0701       {
0702         break;
0703       }
0704       if (direction < 0 && edge_layer <= 7)
0705       {
0706         break;
0707       }
0708 
0709       const unsigned int target_layer =
0710           static_cast<unsigned int>(static_cast<int>(edge_layer) + direction);
0711 
0712       double pred_pad = d->blobs[chain.back()].pad;
0713       double pred_tbin = d->blobs[chain.back()].tbin;
0714       if (direction < 0)
0715       {
0716         pred_pad = d->blobs[chain.front()].pad;
0717         pred_tbin = d->blobs[chain.front()].tbin;
0718       }
0719 
0720       if (chain.size() >= 2)
0721       {
0722         InModuleThreadData::Track tmp;
0723         if (fit_track_from_blobs(d->blobs, chain,
0724                                  d->weight_power, d->adc_weight_floor_frac, tmp))
0725         {
0726           pred_pad = tmp.pad_slope * static_cast<double>(target_layer) + tmp.pad_intercept;
0727           pred_tbin = tmp.tbin_slope * static_cast<double>(target_layer) + tmp.tbin_intercept;
0728         }
0729       }
0730 
0731       const int ibest = find_best_blob_on_layer(d, target_layer, pred_pad, pred_tbin);
0732       if (ibest < 0)
0733       {
0734         break;
0735       }
0736 
0737       d->blobs[ibest].used = 1;
0738       if (direction > 0)
0739       {
0740         chain.push_back(static_cast<unsigned int>(ibest));
0741       }
0742       else
0743       {
0744         chain.insert(chain.begin(), static_cast<unsigned int>(ibest));
0745       }
0746     }
0747   }
0748 
0749   void build_tracks_linear(InModuleThreadData* d)
0750   {
0751     d->tracks.clear();
0752 
0753     std::vector<unsigned int> order;
0754     order.reserve(d->blobs.size());
0755     for (unsigned int i = 0; i < d->blobs.size(); ++i)
0756     {
0757       order.push_back(i);
0758     }
0759     std::sort(order.begin(), order.end(), BlobAdcSort(&d->blobs));
0760 
0761     unsigned int tid = 0;
0762     for (unsigned int seed : order)
0763     {
0764       if (d->blobs[seed].used)
0765       {
0766         continue;
0767       }
0768 
0769       std::vector<unsigned int> chain;
0770       chain.push_back(seed);
0771       d->blobs[seed].used = 1;
0772 
0773       grow_one_direction(d, chain, +1);
0774       grow_one_direction(d, chain, -1);
0775 
0776       if (chain.size() < d->min_track_blobs)
0777       {
0778         for (unsigned int k : chain)
0779         {
0780           d->blobs[k].used = 0;
0781         }
0782         continue;
0783       }
0784 
0785       InModuleThreadData::Track trk;
0786       trk.track_id = tid;
0787       if (make_track_from_blob_chain(d->raw_hits, d->blobs, chain,
0788                                      d->weight_power,
0789                                      d->adc_weight_floor_frac,
0790                                      trk))
0791       {
0792         trk.track_id = tid;
0793         d->tracks.push_back(trk);
0794         ++tid;
0795       }
0796     }
0797   }
0798 
0799   void* ProcessModule(void* arg)
0800   {
0801     InModuleThreadData* d = static_cast<InModuleThreadData*>(arg);
0802     if (!d)
0803     {
0804       return nullptr;
0805     }
0806 
0807     collect_raw_hits(d);
0808     reject_long_pad_noise(d);
0809     build_blobs(d);
0810     build_tracks_linear(d);
0811     connect_track_pieces_in_module(d);
0812 
0813     if (d->verbosity > 1)
0814     {
0815       std::cout << "Tpc_ModuleTrackReco worker: region=" << d->region
0816                 << " sector=" << d->sector << " side=" << d->side
0817                 << " raw_hits=" << d->raw_hits.size()
0818                 << " blobs=" << d->blobs.size()
0819                 << " tracks=" << d->tracks.size() << std::endl;
0820     }
0821 
0822     return nullptr;
0823   }
0824 
0825 }  // anonymous namespace
0826 
0827 // ===================================================================
0828 // Struct constructors
0829 // ===================================================================
0830 InModuleThreadData::LayerHitSet::LayerHitSet()
0831   : layer(0)
0832   , hitsetkey(0)
0833   , hitset(nullptr)
0834 {
0835 }
0836 
0837 InModuleThreadData::RawHit::RawHit()
0838   : layer(0)
0839   , hitsetkey(0)
0840   , hitkey(0)
0841   , pad(0)
0842   , tbin(0)
0843   , adc(0)
0844 {
0845 }
0846 
0847 InModuleThreadData::Blob::Blob()
0848   : layer(0)
0849   , pad(0.0)
0850   , tbin(0.0)
0851   , adc(0.0)
0852   , nhits(0)
0853   , used(0)
0854 {
0855 }
0856 
0857 InModuleThreadData::Track::Track()
0858   : track_id(0)
0859   , first_layer(0)
0860   , last_layer(0)
0861   , nblobs(0)
0862   , nrawhits(0)
0863   , pad_slope(0.0)
0864   , pad_intercept(0.0)
0865   , tbin_slope(0.0)
0866   , tbin_intercept(0.0)
0867 {
0868 }
0869 
0870 InModuleThreadData::InModuleThreadData()
0871   : region(0)
0872   , sector(0)
0873   , side(0)
0874   , module_key(0)
0875   , pedestal(74.4)
0876   , verbosity(0)
0877   , noise_max_consecutive_timebins(10)
0878   , noise_keep_first_timebins(3)
0879   , noise_adc_tolerance(5)
0880   , blob_dt(2)
0881   , blob_dp(2)
0882   , search_dt(6)
0883   , search_dp(6)
0884   , min_track_blobs(4)
0885   , connect_max_layer_gap(8)
0886   , connect_dp(8.0)
0887   , connect_dt(8.0)
0888   , connect_dpad_slope(2.0)
0889   , connect_dtbin_slope(2.0)
0890   , weight_power(0.5)
0891   , adc_weight_floor_frac(0.15)
0892 {
0893 }
0894 
0895 // ===================================================================
0896 // Tpc_ModuleTrackReco
0897 // ===================================================================
0898 Tpc_ModuleTrackReco::Tpc_ModuleTrackReco(const std::string& name,
0899                                          const std::string& filename)
0900   : SubsysReco(name)
0901   , m_outputFileName(filename)
0902   , m_outputFile(nullptr)
0903   , m_tree(nullptr)
0904   , m_hits(nullptr)
0905   , m_tpcModuleTrackContainer(nullptr)
0906   , m_event(0)
0907   , m_maxThreads(72)
0908   , m_pedestal(0.0)
0909   , m_noiseMaxConsecutiveTimebins(10)
0910   , m_noiseKeepFirstTimebins(3)
0911   , m_noiseAdcTolerance(5)
0912   , m_blob_dt(2)
0913   , m_blob_dp(2)
0914   , m_search_dt(6)
0915   , m_search_dp(6)
0916   , m_minTrackBlobs(4)
0917   , m_connectMaxLayerGap(8)
0918   , m_connect_dp(8.0)
0919   , m_connect_dt(8.0)
0920   , m_connect_dpad_slope(2.0)
0921   , m_connect_dtbin_slope(2.0)
0922   , m_tree_event(0)
0923 {
0924 }
0925 
0926 Tpc_ModuleTrackReco::~Tpc_ModuleTrackReco()
0927 {
0928   if (m_outputFile)
0929   {
0930     m_outputFile->Close();
0931     delete m_outputFile;
0932     m_outputFile = nullptr;
0933   }
0934 }
0935 
0936 void Tpc_ModuleTrackReco::setMaxThreads(unsigned int n)
0937 {
0938   m_maxThreads = (n == 0) ? 1 : n;
0939 }
0940 
0941 int Tpc_ModuleTrackReco::Init(PHCompositeNode* /*unused*/)
0942 {
0943   if (Verbosity() <= 0)
0944   {
0945     return Fun4AllReturnCodes::EVENT_OK;
0946   }
0947 
0948   m_outputFile = new TFile(m_outputFileName.c_str(), "RECREATE");
0949   if (!m_outputFile || m_outputFile->IsZombie())
0950   {
0951     std::cerr << Name() << "::Init - cannot create " << m_outputFileName << std::endl;
0952     return Fun4AllReturnCodes::ABORTRUN;
0953   }
0954 
0955   m_tree = new TTree("Tpc_ModuleTrackReco", "TPC in-module pattern recognition");
0956   m_tree->Branch("event", &m_tree_event, "event/I");
0957 
0958   m_tree->Branch("track_id", &m_tree_track_id);
0959   m_tree->Branch("region", &m_tree_region);
0960   m_tree->Branch("sector", &m_tree_sector);
0961   m_tree->Branch("side", &m_tree_side);
0962   m_tree->Branch("nblobs", &m_tree_nblobs);
0963   m_tree->Branch("nrawhits", &m_tree_nrawhits);
0964   m_tree->Branch("first_layer", &m_tree_first_layer);
0965   m_tree->Branch("last_layer", &m_tree_last_layer);
0966 
0967   // Per-hit branches: TrkrHitSetContainer keys only
0968   m_tree->Branch("hit_event", &m_tree_hit_event);
0969   m_tree->Branch("hit_track_id", &m_tree_hit_track_id);
0970   m_tree->Branch("hit_region", &m_tree_hit_region);
0971   m_tree->Branch("hit_sector", &m_tree_hit_sector);
0972   m_tree->Branch("hit_side", &m_tree_hit_side);
0973   m_tree->Branch("hit_layer", &m_tree_hit_layer);
0974   m_tree->Branch("hit_hitsetkey", &m_tree_hit_hitsetkey);
0975   m_tree->Branch("hit_hitkey", &m_tree_hit_hitkey);
0976 
0977   std::cout << Name() << "::Init - output file " << m_outputFileName << " created" << std::endl;
0978 
0979   return Fun4AllReturnCodes::EVENT_OK;
0980 }
0981 
0982 int Tpc_ModuleTrackReco::InitRun(PHCompositeNode* topNode)
0983 {
0984   if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0985   {
0986     return Fun4AllReturnCodes::ABORTRUN;
0987   }
0988   if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0989   {
0990     return Fun4AllReturnCodes::ABORTRUN;
0991   }
0992 
0993   m_event = 0;
0994   return Fun4AllReturnCodes::EVENT_OK;
0995 }
0996 
0997 int Tpc_ModuleTrackReco::End(PHCompositeNode* /*unused*/)
0998 {
0999   if (m_outputFile)
1000   {
1001     m_outputFile->cd();
1002     if (m_tree)
1003     {
1004       m_tree->Write();
1005     }
1006     m_outputFile->Close();
1007     delete m_outputFile;
1008     m_outputFile = nullptr;
1009   }
1010   return Fun4AllReturnCodes::EVENT_OK;
1011 }
1012 
1013 int Tpc_ModuleTrackReco::getNodes(PHCompositeNode* topNode)
1014 {
1015   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1016   if (!m_hits)
1017   {
1018     std::cerr << Name() << "::getNodes - missing TRKR_HITSET" << std::endl;
1019     return Fun4AllReturnCodes::ABORTRUN;
1020   }
1021 
1022   return Fun4AllReturnCodes::EVENT_OK;
1023 }
1024 
1025 int Tpc_ModuleTrackReco::createNodes(PHCompositeNode* topNode)
1026 {
1027   PHNodeIterator iter(topNode);
1028 
1029   PHCompositeNode* dstNode =
1030       dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
1031 
1032   if (!dstNode)
1033   {
1034     dstNode = new PHCompositeNode("DST");
1035     topNode->addNode(dstNode);
1036   }
1037 
1038   m_tpcModuleTrackContainer =
1039       findNode::getClass<Tpc_ModuleTrackContainer>(topNode, "TPC_MODULETRACKS");
1040 
1041   if (!m_tpcModuleTrackContainer)
1042   {
1043     m_tpcModuleTrackContainer = new Tpc_ModuleTrackContainerv1();
1044 
1045     PHIODataNode<PHObject>* node =
1046         new PHIODataNode<PHObject>(m_tpcModuleTrackContainer,
1047                                    "TPC_MODULETRACKS", "PHObject");
1048     dstNode->addNode(node);
1049 
1050     std::cout << Name() << "::createNodes - created TPC_MODULETRACKS node" << std::endl;
1051   }
1052 
1053   return Fun4AllReturnCodes::EVENT_OK;
1054 }
1055 
1056 void Tpc_ModuleTrackReco::reset_tree_vars()
1057 {
1058   m_tree_event = m_event;
1059 
1060   m_tree_track_id.clear();
1061   m_tree_region.clear();
1062   m_tree_sector.clear();
1063   m_tree_side.clear();
1064   m_tree_nblobs.clear();
1065   m_tree_nrawhits.clear();
1066   m_tree_first_layer.clear();
1067   m_tree_last_layer.clear();
1068 
1069   m_tree_hit_event.clear();
1070   m_tree_hit_track_id.clear();
1071   m_tree_hit_region.clear();
1072   m_tree_hit_sector.clear();
1073   m_tree_hit_side.clear();
1074   m_tree_hit_layer.clear();
1075   m_tree_hit_hitsetkey.clear();
1076   m_tree_hit_hitkey.clear();
1077 }
1078 
1079 int Tpc_ModuleTrackReco::process_event(PHCompositeNode* /*unused*/)
1080 {
1081   reset_tree_vars();
1082 
1083   if (m_tpcModuleTrackContainer)
1084   {
1085     m_tpcModuleTrackContainer->Reset();
1086   }
1087 
1088   std::vector<InModuleThreadData> tdata;
1089   tdata.reserve(72);
1090 
1091   for (unsigned int side = 0; side < 2; ++side)
1092   {
1093     for (unsigned int sector = 0; sector < 12; ++sector)
1094     {
1095       for (unsigned int region = 0; region < 3; ++region)
1096       {
1097         InModuleThreadData td;
1098         td.region = region;
1099         td.sector = sector;
1100         td.side = static_cast<int>(side);
1101         td.module_key = TpcDefs::genModuleHitSetKey(static_cast<uint8_t>(region),
1102                                                     static_cast<uint8_t>(sector),
1103                                                     static_cast<uint8_t>(side));
1104         td.pedestal = m_pedestal;
1105         td.verbosity = Verbosity();
1106         td.noise_max_consecutive_timebins = m_noiseMaxConsecutiveTimebins;
1107         td.noise_keep_first_timebins = m_noiseKeepFirstTimebins;
1108         td.noise_adc_tolerance = m_noiseAdcTolerance;
1109         td.blob_dt = m_blob_dt;
1110         td.blob_dp = m_blob_dp;
1111         td.search_dt = m_search_dt;
1112         td.search_dp = m_search_dp;
1113         td.min_track_blobs = m_minTrackBlobs;
1114         td.connect_max_layer_gap = m_connectMaxLayerGap;
1115         td.connect_dp = m_connect_dp;
1116         td.connect_dt = m_connect_dt;
1117         td.connect_dpad_slope = m_connect_dpad_slope;
1118         td.connect_dtbin_slope = m_connect_dtbin_slope;
1119 
1120         for (unsigned int l = 0; l < 16; ++l)
1121         {
1122           const unsigned int layer = region * 16 + l + 7;
1123           const TrkrDefs::hitsetkey hitset_key = TpcDefs::genHitSetKey(layer, sector, side);
1124           TrkrHitSet* hitset = m_hits->findHitSet(hitset_key);
1125           if (!hitset)
1126           {
1127             continue;
1128           }
1129 
1130           InModuleThreadData::LayerHitSet lhs;
1131           lhs.layer = layer;
1132           lhs.hitsetkey = hitset_key;
1133           lhs.hitset = hitset;
1134           td.layer_hitsets.push_back(lhs);
1135         }
1136 
1137         if (!td.layer_hitsets.empty())
1138         {
1139           tdata.push_back(td);
1140         }
1141       }
1142     }
1143   }
1144   if(Verbosity() > 1)
1145   {
1146   std::cout << Name() << "::process_event - event " << m_event
1147             << " has " << tdata.size() << " non-empty modules" << std::endl;
1148   }
1149   const unsigned int maxLive = std::max(1U,
1150                                         std::min(m_maxThreads, static_cast<unsigned int>(tdata.size())));
1151 
1152   for (unsigned int start = 0; start < static_cast<unsigned int>(tdata.size()); start += maxLive)
1153   {
1154     const unsigned int end = std::min(start + maxLive,
1155                                       static_cast<unsigned int>(tdata.size()));
1156     const unsigned int nLive = end - start;
1157 
1158     std::vector<pthread_t> threads(nLive);
1159     std::vector<int> thread_ok(nLive, 0);
1160 
1161     for (unsigned int i = 0; i < nLive; ++i)
1162     {
1163       const unsigned int idx = start + i;
1164       const int rc = pthread_create(&threads[i], nullptr, ProcessModule,
1165                                     static_cast<void*>(&tdata[idx]));
1166       if (rc != 0)
1167       {
1168         std::cerr << Name() << "::process_event - pthread_create failed for"
1169                   << " region=" << tdata[idx].region
1170                   << " sector=" << tdata[idx].sector
1171                   << " side=" << tdata[idx].side << std::endl;
1172       }
1173       else
1174       {
1175         thread_ok[i] = 1;
1176       }
1177     }
1178 
1179     for (unsigned int i = 0; i < nLive; ++i)
1180     {
1181       if (thread_ok[i])
1182       {
1183         pthread_join(threads[i], nullptr);
1184       }
1185     }
1186   }
1187 
1188   // Harvest results from all modules
1189   for (const auto& td : tdata)
1190   {
1191     for (unsigned int it = 0; it < td.tracks.size(); ++it)
1192     {
1193       const InModuleThreadData::Track& tr = td.tracks[it];
1194 
1195       const unsigned int global_track_id =
1196           m_tpcModuleTrackContainer ? m_tpcModuleTrackContainer->size() : static_cast<unsigned int>(m_tree_track_id.size());
1197 
1198       Tpc_ModuleTrackv1* outTrack = new Tpc_ModuleTrackv1();
1199 
1200       outTrack->set_event(static_cast<unsigned int>(m_event));
1201       outTrack->set_track_id(global_track_id);
1202       outTrack->set_region(td.region);
1203       outTrack->set_sector(td.sector);
1204       outTrack->set_side(td.side);
1205       outTrack->set_nblobs(tr.nblobs);
1206       outTrack->set_nrawhits(tr.nrawhits);
1207       outTrack->set_first_layer(tr.first_layer);
1208       outTrack->set_last_layer(tr.last_layer);
1209 
1210       // TTree track-level fill
1211       m_tree_track_id.push_back(global_track_id);
1212       m_tree_region.push_back(td.region);
1213       m_tree_sector.push_back(td.sector);
1214       m_tree_side.push_back(td.side);
1215       m_tree_nblobs.push_back(tr.nblobs);
1216       m_tree_nrawhits.push_back(tr.nrawhits);
1217       m_tree_first_layer.push_back(tr.first_layer);
1218       m_tree_last_layer.push_back(tr.last_layer);
1219 
1220       // Hit references: (hitsetkey, hitkey) only — no data copy
1221       for (unsigned int raw_hit_indice : tr.raw_hit_indices)
1222       {
1223         const InModuleThreadData::RawHit& rh = td.raw_hits[raw_hit_indice];
1224 
1225         outTrack->add_hit_index(rh.hitsetkey, rh.hitkey);
1226 
1227         m_tree_hit_event.push_back(static_cast<unsigned int>(m_event));
1228         m_tree_hit_track_id.push_back(global_track_id);
1229         m_tree_hit_region.push_back(td.region);
1230         m_tree_hit_sector.push_back(td.sector);
1231         m_tree_hit_side.push_back(td.side);
1232         m_tree_hit_layer.push_back(rh.layer);
1233         m_tree_hit_hitsetkey.push_back(static_cast<unsigned long long>(rh.hitsetkey));
1234         m_tree_hit_hitkey.push_back(static_cast<unsigned long long>(rh.hitkey));
1235       }
1236 
1237       if (m_tpcModuleTrackContainer)
1238       {
1239         m_tpcModuleTrackContainer->add_track(outTrack);
1240       }
1241       else
1242       {
1243         delete outTrack;
1244       }
1245     }
1246   }
1247 
1248   if (m_tree)
1249   {
1250     m_tree->Fill();
1251   }
1252 
1253   if (Verbosity() > 0)
1254   {
1255     std::cout << Name() << "::process_event - event " << m_event
1256               << " tracks=" << m_tree_track_id.size()
1257               << " track-raw-hits=" << m_tree_hit_track_id.size() << std::endl;
1258   }
1259   if (m_tpcModuleTrackContainer && Verbosity() > 0)
1260   {
1261     m_tpcModuleTrackContainer->identify();
1262   }
1263   if (m_tpcModuleTrackContainer && Verbosity() > 1)
1264   {
1265     for (unsigned int i = 0; i < m_tpcModuleTrackContainer->size(); ++i)
1266     {
1267       const Tpc_ModuleTrack* trk = m_tpcModuleTrackContainer->get_track(i);
1268       if (trk)
1269       {
1270         trk->identify();
1271       }
1272     }
1273   }
1274 
1275   ++m_event;
1276   return Fun4AllReturnCodes::EVENT_OK;
1277 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!