Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-07-16 08:14:51 

0001 #include "PHTruthTrackFitter.h"
0002 
0003 #include <trackbase/ActsGeometry.h>
0004 #include <trackbase/TrkrCluster.h>
0005 #include <trackbase/TrkrClusterContainer.h>
0006 #include <trackbase/TrkrClusterHitAssoc.h>
0007 #include <trackbase/TrkrHitTruthAssoc.h>
0008 
0009 #include <trackbase_historic/SvtxTrack.h>
0010 #include <trackbase_historic/SvtxTrackMap.h>
0011 #include <trackbase_historic/SvtxTrackMap_v2.h>
0012 #include <trackbase_historic/SvtxTrackState_v3.h>
0013 #include <trackbase_historic/SvtxTrack_v4.h>
0014 #include <trackbase_historic/TrackSeed.h>
0015 #include <trackbase_historic/TrackSeedContainer.h>
0016 
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018 
0019 #include <g4main/PHG4Hit.h>
0020 #include <g4main/PHG4HitContainer.h>
0021 #include <g4main/PHG4Particle.h>
0022 #include <g4main/PHG4TruthInfoContainer.h>
0023 #include <g4main/PHG4VtxPoint.h>
0024 
0025 #include <phool/PHCompositeNode.h>
0026 #include <phool/PHIODataNode.h>
0027 #include <phool/PHNode.h>
0028 #include <phool/PHNodeIterator.h>
0029 #include <phool/PHObject.h>
0030 #include <phool/getClass.h>
0031 #include <phool/phool.h>
0032 
0033 #include <TDatabasePDG.h>
0034 #include <TParticlePDG.h>
0035 
0036 #include <algorithm>
0037 #include <cmath>
0038 #include <iostream>
0039 #include <limits>
0040 #include <set>
0041 #include <utility>
0042 #include <vector>
0043 
0044 namespace
0045 {
0046   template <class T>
0047   constexpr T square(const T& x)
0048   {
0049     return x * x;
0050   }
0051 
0052   bool is_finite(float value)
0053   {
0054     return std::isfinite(value);
0055   }
0056 
0057   float average_or(float a, float b, float fallback)
0058   {
0059     const bool aok = is_finite(a);
0060     const bool bok = is_finite(b);
0061 
0062     if (aok && bok)
0063     {
0064       return 0.5 * (a + b);
0065     }
0066     if (aok)
0067     {
0068       return a;
0069     }
0070     if (bok)
0071     {
0072       return b;
0073     }
0074 
0075     return fallback;
0076   }
0077 
0078   bool valid_track_id(unsigned int trackid)
0079   {
0080     return trackid != std::numeric_limits<unsigned int>::max();
0081   }
0082 
0083   class InterpolationData
0084   {
0085    public:
0086     InterpolationData(double x, double y, double z, double px, double py, double pz, double weight)
0087       : m_x(x)
0088       , m_y(y)
0089       , m_z(z)
0090       , m_px(px)
0091       , m_py(py)
0092       , m_pz(pz)
0093       , m_weight(weight)
0094     {
0095     }
0096 
0097     double r() const
0098     {
0099       return std::sqrt(square(m_x) + square(m_y));
0100     }
0101 
0102     double x() const { return m_x; }
0103     double y() const { return m_y; }
0104     double z() const { return m_z; }
0105     double px() const { return m_px; }
0106     double py() const { return m_py; }
0107     double pz() const { return m_pz; }
0108     double weight() const { return m_weight; }
0109 
0110    private:
0111     double m_x = 0;
0112     double m_y = 0;
0113     double m_z = 0;
0114     double m_px = 0;
0115     double m_py = 0;
0116     double m_pz = 0;
0117     double m_weight = 1;
0118   };
0119 
0120   template <double (InterpolationData::*accessor)() const>
0121   double interpolate_r(const std::vector<InterpolationData>& hits, double r_extrap, double fallback)
0122   {
0123     double sw = 0;
0124     double swr = 0;
0125     double swr2 = 0;
0126     double swq = 0;
0127     double swrq = 0;
0128 
0129     for (const auto& hit : hits)
0130     {
0131       const auto q = (hit.*accessor)();
0132       const auto r = hit.r();
0133       const auto weight = hit.weight();
0134       if (!std::isfinite(q) || !std::isfinite(r) || !std::isfinite(weight) || weight <= 0)
0135       {
0136         continue;
0137       }
0138 
0139       sw += weight;
0140       swr += weight * r;
0141       swr2 += weight * square(r);
0142       swq += weight * q;
0143       swrq += weight * r * q;
0144     }
0145 
0146     /*
0147      * Fit q(r) = a*r + b with weighted least squares, where q is one of
0148      * x/y/z/px/py/pz. The sums above form the normal equations:
0149      *
0150      *   a*swr2 + b*swr = swrq
0151      *   a*swr  + b*sw  = swq
0152      *
0153      * alpha and beta are the Cramer's-rule numerators for the slope and
0154      * intercept. Keeping the final division common is the same as returning
0155      * slope*r_extrap + intercept, but avoids one extra division.
0156      */
0157     const auto denom = sw * swr2 - square(swr);
0158     const auto scale = std::max(std::abs(sw * swr2), square(swr));
0159     if (scale <= 0 || std::abs(denom) <= std::numeric_limits<double>::epsilon() * scale)
0160     {
0161       return fallback;
0162     }
0163 
0164     const auto alpha = sw * swrq - swr * swq;
0165     const auto beta = swr2 * swq - swr * swrq;
0166     const auto value = (alpha * r_extrap + beta) / denom;
0167     return std::isfinite(value) ? value : fallback;
0168   }
0169 }  // namespace
0170 
0171 PHTruthTrackFitter::PHTruthTrackFitter(const std::string& name)
0172   : SubsysReco(name)
0173 {
0174 }
0175 
0176 int PHTruthTrackFitter::InitRun(PHCompositeNode* topNode)
0177 {
0178   if (Verbosity() > 0)
0179   {
0180     std::cout << "PHTruthTrackFitter::InitRun - output track map: " << m_trackMapName << std::endl;
0181   }
0182 
0183   if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0184   {
0185     return Fun4AllReturnCodes::ABORTEVENT;
0186   }
0187 
0188   if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0189   {
0190     return Fun4AllReturnCodes::ABORTEVENT;
0191   }
0192 
0193   return Fun4AllReturnCodes::EVENT_OK;
0194 }
0195 
0196 int PHTruthTrackFitter::process_event(PHCompositeNode* /*topNode*/)
0197 {
0198   m_trackMap->Reset();
0199 
0200   unsigned int skipped_tracks = 0;
0201   for (auto* seed : *m_seedMap)
0202   {
0203     if (!seed)
0204     {
0205       continue;
0206     }
0207 
0208     auto* tpc_seed = getSeed(m_tpcSeeds, seed->get_tpc_seed_index());
0209     auto* silicon_seed = getSeed(m_siliconSeeds, seed->get_silicon_seed_index());
0210 
0211     if (!tpc_seed && !silicon_seed)
0212     {
0213       ++skipped_tracks;
0214       continue;
0215     }
0216 
0217     const auto truth_track_id = getTruthTrackId(seed, tpc_seed, silicon_seed);
0218     if (!valid_track_id(truth_track_id))
0219     {
0220       if (Verbosity() > 1)
0221       {
0222         std::cout << "PHTruthTrackFitter::process_event - could not determine truth id for seed" << std::endl;
0223       }
0224       ++skipped_tracks;
0225       continue;
0226     }
0227 
0228     auto* g4particle = m_g4TruthInfo->GetParticle(truth_track_id);
0229     if (!g4particle)
0230     {
0231       if (Verbosity() > 1)
0232       {
0233         std::cout << "PHTruthTrackFitter::process_event - no PHG4Particle for track id "
0234                   << truth_track_id << std::endl;
0235       }
0236       ++skipped_tracks;
0237       continue;
0238     }
0239 
0240     const auto* g4vertex = m_g4TruthInfo->GetVtx(g4particle->get_vtx_id());
0241     if (!g4vertex)
0242     {
0243       if (Verbosity() > 1)
0244       {
0245         std::cout << "PHTruthTrackFitter::process_event - no PHG4VtxPoint for track id "
0246                   << truth_track_id << std::endl;
0247       }
0248       ++skipped_tracks;
0249       continue;
0250     }
0251 
0252     SvtxTrack_v4 track;
0253     track.set_tpc_seed(tpc_seed);
0254     track.set_silicon_seed(silicon_seed);
0255     track.set_crossing(getCrossing(tpc_seed, silicon_seed));
0256     track.set_vertex_id(g4particle->get_vtx_id());
0257     track.set_charge(getCharge(g4particle, tpc_seed, silicon_seed));
0258     track.set_chisq(0);
0259 
0260     track.set_x(g4vertex->get_x());
0261     track.set_y(g4vertex->get_y());
0262     track.set_z(g4vertex->get_z());
0263     track.set_px(g4particle->get_px());
0264     track.set_py(g4particle->get_py());
0265     track.set_pz(g4particle->get_pz());
0266 
0267     for (int i = 0; i < 6; ++i)
0268     {
0269       for (int j = i; j < 6; ++j)
0270       {
0271         track.set_error(i, j, 0);
0272       }
0273     }
0274     track.set_error(0, 0, square(m_positionError));
0275     track.set_error(1, 1, square(m_positionError));
0276     track.set_error(2, 2, square(m_zError));
0277 
0278     unsigned int state_index = 1;
0279     for (const auto* track_seed : {silicon_seed, tpc_seed})
0280     {
0281       if (!track_seed)
0282       {
0283         continue;
0284       }
0285 
0286       for (auto iter = track_seed->begin_cluster_keys(); iter != track_seed->end_cluster_keys(); ++iter)
0287       {
0288         if (addStateFromCluster(&track, *iter, truth_track_id, g4particle, g4vertex, state_index))
0289         {
0290           ++state_index;
0291         }
0292       }
0293     }
0294 
0295     if (track.size_states() <= 1)
0296     {
0297       if (Verbosity() > 1)
0298       {
0299         std::cout << "PHTruthTrackFitter::process_event - no truth states for track id "
0300                   << truth_track_id << std::endl;
0301       }
0302       ++skipped_tracks;
0303       continue;
0304     }
0305 
0306     track.set_ndf(std::max<int>(0, 2 * static_cast<int>(track.size_states()) - 5));
0307 
0308     const unsigned int track_id = m_trackMap->size();
0309     track.set_id(track_id);
0310     m_trackMap->insertWithKey(&track, track_id);
0311   }
0312 
0313   if (Verbosity() > 0)
0314   {
0315     std::cout << "PHTruthTrackFitter::process_event - built " << m_trackMap->size()
0316               << " truth tracks, skipped " << skipped_tracks << std::endl;
0317   }
0318 
0319   return Fun4AllReturnCodes::EVENT_OK;
0320 }
0321 
0322 int PHTruthTrackFitter::End(PHCompositeNode* /*topNode*/)
0323 {
0324   return Fun4AllReturnCodes::EVENT_OK;
0325 }
0326 
0327 int PHTruthTrackFitter::createNodes(PHCompositeNode* topNode)
0328 {
0329   PHNodeIterator iter(topNode);
0330 
0331   auto* dst_node = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0332   if (!dst_node)
0333   {
0334     std::cerr << PHWHERE << "DST node is missing" << std::endl;
0335     return Fun4AllReturnCodes::ABORTEVENT;
0336   }
0337 
0338   PHNodeIterator dst_iter(dst_node);
0339   auto* svtx_node = dynamic_cast<PHCompositeNode*>(dst_iter.findFirst("PHCompositeNode", "SVTX"));
0340   if (!svtx_node)
0341   {
0342     svtx_node = new PHCompositeNode("SVTX");
0343     dst_node->addNode(svtx_node);
0344   }
0345 
0346   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0347   if (!m_trackMap)
0348   {
0349     m_trackMap = new SvtxTrackMap_v2;
0350     auto* track_node = new PHIODataNode<PHObject>(m_trackMap, m_trackMapName, "PHObject");
0351     svtx_node->addNode(track_node);
0352   }
0353 
0354   return Fun4AllReturnCodes::EVENT_OK;
0355 }
0356 
0357 int PHTruthTrackFitter::getNodes(PHCompositeNode* topNode)
0358 {
0359   m_seedMap = findNode::getClass<TrackSeedContainer>(topNode, m_svtxSeedMapName);
0360   if (!m_seedMap)
0361   {
0362     std::cout << PHWHERE << "No " << m_svtxSeedMapName << " on node tree. Bailing" << std::endl;
0363     return Fun4AllReturnCodes::ABORTEVENT;
0364   }
0365 
0366   m_tpcSeeds = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0367   if (!m_tpcSeeds)
0368   {
0369     std::cout << PHWHERE << "No TpcTrackSeedContainer on node tree. Bailing" << std::endl;
0370     return Fun4AllReturnCodes::ABORTEVENT;
0371   }
0372 
0373   m_siliconSeeds = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0374   if (!m_siliconSeeds)
0375   {
0376     std::cout << PHWHERE << "No SiliconTrackSeedContainer on node tree. Bailing" << std::endl;
0377     return Fun4AllReturnCodes::ABORTEVENT;
0378   }
0379 
0380   m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterMapName);
0381   if (!m_clusterMap)
0382   {
0383     std::cout << PHWHERE << "No " << m_clusterMapName << " on node tree. Bailing" << std::endl;
0384     return Fun4AllReturnCodes::ABORTEVENT;
0385   }
0386 
0387   m_clusterHitMap = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0388   if (!m_clusterHitMap)
0389   {
0390     std::cout << PHWHERE << "No TRKR_CLUSTERHITASSOC on node tree. Bailing" << std::endl;
0391     return Fun4AllReturnCodes::ABORTEVENT;
0392   }
0393 
0394   m_hitTruthAssoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0395   if (!m_hitTruthAssoc)
0396   {
0397     std::cout << PHWHERE << "No TRKR_HITTRUTHASSOC on node tree. Bailing" << std::endl;
0398     return Fun4AllReturnCodes::ABORTEVENT;
0399   }
0400 
0401   m_g4TruthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0402   if (!m_g4TruthInfo)
0403   {
0404     std::cout << PHWHERE << "No G4TruthInfo on node tree. Bailing" << std::endl;
0405     return Fun4AllReturnCodes::ABORTEVENT;
0406   }
0407 
0408   m_g4HitsTpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
0409   m_g4HitsIntt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
0410   m_g4HitsMvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
0411   m_g4HitsMicromegas = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
0412 
0413   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0414   if (m_extrapolateToClusterRadius && !m_tGeometry)
0415   {
0416     std::cout << PHWHERE << "No ActsGeometry on node tree. Bailing" << std::endl;
0417     return Fun4AllReturnCodes::ABORTEVENT;
0418   }
0419 
0420   return Fun4AllReturnCodes::EVENT_OK;
0421 }
0422 
0423 TrackSeed* PHTruthTrackFitter::getSeed(TrackSeedContainer* container, unsigned int index) const
0424 {
0425   if (!container || index >= container->size())
0426   {
0427     return nullptr;
0428   }
0429 
0430   return container->get(index);
0431 }
0432 
0433 unsigned int PHTruthTrackFitter::getTruthTrackId(const TrackSeed* svtxSeed,
0434                                                  const TrackSeed* tpcSeed,
0435                                                  const TrackSeed* siliconSeed) const
0436 {
0437   auto truth_track_id = m_invalidTruthTrackId;
0438   for (const auto* seed : {svtxSeed, tpcSeed, siliconSeed})
0439   {
0440     if (!seed)
0441     {
0442       continue;
0443     }
0444 
0445     const auto seed_truth_track_id = seed->get_truth_track_id();
0446     if (!valid_track_id(seed_truth_track_id))
0447     {
0448       continue;
0449     }
0450 
0451     if (valid_track_id(truth_track_id) && seed_truth_track_id != truth_track_id)
0452     {
0453       if (Verbosity() > 0)
0454       {
0455         std::cout << "PHTruthTrackFitter::getTruthTrackId - inconsistent seed truth ids "
0456                   << truth_track_id << " and " << seed_truth_track_id << std::endl;
0457       }
0458       return m_invalidTruthTrackId;
0459     }
0460 
0461     truth_track_id = seed_truth_track_id;
0462   }
0463 
0464   return truth_track_id;
0465 }
0466 
0467 std::vector<const PHG4Hit*> PHTruthTrackFitter::getTruthHits(TrkrDefs::cluskey cluskey) const
0468 {
0469   std::vector<const PHG4Hit*> truth_hits;
0470   if (!m_clusterHitMap || !m_hitTruthAssoc)
0471   {
0472     return truth_hits;
0473   }
0474 
0475   const auto hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluskey);
0476   const auto trkrid = TrkrDefs::getTrkrId(hitsetkey);
0477   const auto hitrange = m_clusterHitMap->getHits(cluskey);
0478 
0479   std::set<PHG4HitDefs::keytype> used_g4hits;
0480   for (auto clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
0481   {
0482     const auto hitkey = clushititer->second;
0483 
0484     TrkrHitTruthAssoc::MMap temp_map;
0485     m_hitTruthAssoc->getG4Hits(hitsetkey, hitkey, temp_map);
0486 
0487     for (const auto& hit_truth_iter : temp_map)
0488     {
0489       const auto g4hitkey = hit_truth_iter.second.second;
0490       if (!used_g4hits.insert(g4hitkey).second)
0491       {
0492         continue;
0493       }
0494 
0495       const auto* g4hit = getG4Hit(trkrid, g4hitkey);
0496       if (g4hit)
0497       {
0498         truth_hits.push_back(g4hit);
0499       }
0500     }
0501   }
0502 
0503   return truth_hits;
0504 }
0505 
0506 const PHG4Hit* PHTruthTrackFitter::getG4Hit(unsigned int trkrid, PHG4HitDefs::keytype g4hitkey) const
0507 {
0508   PHG4HitContainer* container = nullptr;
0509   switch (trkrid)
0510   {
0511   case TrkrDefs::tpcId:
0512     container = m_g4HitsTpc;
0513     break;
0514   case TrkrDefs::inttId:
0515     container = m_g4HitsIntt;
0516     break;
0517   case TrkrDefs::mvtxId:
0518     container = m_g4HitsMvtx;
0519     break;
0520   case TrkrDefs::micromegasId:
0521     container = m_g4HitsMicromegas;
0522     break;
0523   default:
0524     break;
0525   }
0526 
0527   return container ? container->findHit(g4hitkey) : nullptr;
0528 }
0529 
0530 bool PHTruthTrackFitter::addStateFromCluster(SvtxTrack* track,
0531                                              TrkrDefs::cluskey cluskey,
0532                                              unsigned int truthTrackId,
0533                                              const PHG4Particle* particle,
0534                                              const PHG4VtxPoint* vertex,
0535                                              unsigned int stateIndex) const
0536 {
0537   auto* cluster = m_clusterMap->findCluster(cluskey);
0538   if (!cluster)
0539   {
0540     return false;
0541   }
0542 
0543   std::vector<InterpolationData> interpolation_hits;
0544   double weight_sum = 0;
0545   double x = 0;
0546   double y = 0;
0547   double z = 0;
0548   double px = 0;
0549   double py = 0;
0550   double pz = 0;
0551   double local_x = 0;
0552   double local_y = 0;
0553 
0554   for (const auto* g4hit : getTruthHits(cluskey))
0555   {
0556     if (!g4hit || g4hit->get_trkid() != static_cast<int>(truthTrackId))
0557     {
0558       continue;
0559     }
0560 
0561     const auto hit_x = average_or(g4hit->get_x(0), g4hit->get_x(1), std::numeric_limits<float>::quiet_NaN());
0562     const auto hit_y = average_or(g4hit->get_y(0), g4hit->get_y(1), std::numeric_limits<float>::quiet_NaN());
0563     const auto hit_z = average_or(g4hit->get_z(0), g4hit->get_z(1), std::numeric_limits<float>::quiet_NaN());
0564     if (!is_finite(hit_x) || !is_finite(hit_y) || !is_finite(hit_z))
0565     {
0566       continue;
0567     }
0568 
0569     const auto hit_px = average_or(g4hit->get_px(0), g4hit->get_px(1), particle->get_px());
0570     const auto hit_py = average_or(g4hit->get_py(0), g4hit->get_py(1), particle->get_py());
0571     const auto hit_pz = average_or(g4hit->get_pz(0), g4hit->get_pz(1), particle->get_pz());
0572     const auto hit_local_x = average_or(g4hit->get_local_x(0), g4hit->get_local_x(1), 0);
0573     const auto hit_local_y = average_or(g4hit->get_local_y(0), g4hit->get_local_y(1), 0);
0574 
0575     double weight = g4hit->get_edep();
0576     if (!std::isfinite(weight) || weight <= 0)
0577     {
0578       weight = 1;
0579     }
0580 
0581     for (int endpoint = 0; endpoint < 2; ++endpoint)
0582     {
0583       const auto endpoint_x = g4hit->get_x(endpoint);
0584       const auto endpoint_y = g4hit->get_y(endpoint);
0585       const auto endpoint_z = g4hit->get_z(endpoint);
0586       if (!is_finite(endpoint_x) || !is_finite(endpoint_y) || !is_finite(endpoint_z))
0587       {
0588         continue;
0589       }
0590 
0591       const auto endpoint_px = g4hit->get_px(endpoint);
0592       const auto endpoint_py = g4hit->get_py(endpoint);
0593       const auto endpoint_pz = g4hit->get_pz(endpoint);
0594 
0595       interpolation_hits.emplace_back(endpoint_x,
0596                                       endpoint_y,
0597                                       endpoint_z,
0598                                       is_finite(endpoint_px) ? endpoint_px : particle->get_px(),
0599                                       is_finite(endpoint_py) ? endpoint_py : particle->get_py(),
0600                                       is_finite(endpoint_pz) ? endpoint_pz : particle->get_pz(),
0601                                       weight);
0602     }
0603 
0604     weight_sum += weight;
0605     x += weight * hit_x;
0606     y += weight * hit_y;
0607     z += weight * hit_z;
0608     px += weight * hit_px;
0609     py += weight * hit_py;
0610     pz += weight * hit_pz;
0611     local_x += weight * hit_local_x;
0612     local_y += weight * hit_local_y;
0613   }
0614 
0615   if (weight_sum <= 0)
0616   {
0617     return false;
0618   }
0619 
0620   x /= weight_sum;
0621   y /= weight_sum;
0622   z /= weight_sum;
0623   px /= weight_sum;
0624   py /= weight_sum;
0625   pz /= weight_sum;
0626   local_x /= weight_sum;
0627   local_y /= weight_sum;
0628 
0629   if (m_extrapolateToClusterRadius && !interpolation_hits.empty())
0630   {
0631     const auto cluster_radius = getClusterRadius(cluskey, cluster);
0632     if (std::isfinite(cluster_radius) && cluster_radius > 0)
0633     {
0634       x = interpolate_r<&InterpolationData::x>(interpolation_hits, cluster_radius, x);
0635       y = interpolate_r<&InterpolationData::y>(interpolation_hits, cluster_radius, y);
0636       z = interpolate_r<&InterpolationData::z>(interpolation_hits, cluster_radius, z);
0637       px = interpolate_r<&InterpolationData::px>(interpolation_hits, cluster_radius, px);
0638       py = interpolate_r<&InterpolationData::py>(interpolation_hits, cluster_radius, py);
0639       pz = interpolate_r<&InterpolationData::pz>(interpolation_hits, cluster_radius, pz);
0640     }
0641     else if (Verbosity() > 1)
0642     {
0643       std::cout << "PHTruthTrackFitter::addStateFromCluster - invalid cluster radius for cluster "
0644                 << cluskey << ", using truth hit average" << std::endl;
0645     }
0646   }
0647 
0648   float pathlength = getPathLength(vertex, x, y, z, stateIndex);
0649   while (track->count_states(pathlength) != 0)
0650   {
0651     pathlength += 1.e-3;
0652   }
0653 
0654   SvtxTrackState_v3 state(pathlength);
0655   state.set_name("PHTruthTrackFitter");
0656   state.set_cluskey(cluskey);
0657   state.set_x(x);
0658   state.set_y(y);
0659   state.set_z(z);
0660   state.set_px(px);
0661   state.set_py(py);
0662   state.set_pz(pz);
0663   state.set_localX(local_x);
0664   state.set_localY(local_y);
0665 
0666   for (int i = 0; i < 6; ++i)
0667   {
0668     for (int j = i; j < 6; ++j)
0669     {
0670       state.set_error(i, j, 0);
0671     }
0672   }
0673   state.set_error(0, 0, square(m_positionError));
0674   state.set_error(1, 1, square(m_positionError));
0675   state.set_error(2, 2, square(m_zError));
0676 
0677   track->insert_state(&state);
0678   return true;
0679 }
0680 
0681 float PHTruthTrackFitter::getClusterRadius(TrkrDefs::cluskey cluskey, TrkrCluster* cluster) const
0682 {
0683   if (!m_tGeometry || !cluster)
0684   {
0685     return std::numeric_limits<float>::quiet_NaN();
0686   }
0687 
0688   const auto global = m_tGeometry->getGlobalPosition(cluskey, cluster);
0689   const auto radius = std::sqrt(square(global.x()) + square(global.y()));
0690   return std::isfinite(radius) ? radius : std::numeric_limits<float>::quiet_NaN();
0691 }
0692 
0693 float PHTruthTrackFitter::getPathLength(const PHG4VtxPoint* vertex,
0694                                         float x, float y, float z,
0695                                         unsigned int stateIndex) const
0696 {
0697   if (vertex)
0698   {
0699     const auto dx = x - vertex->get_x();
0700     const auto dy = y - vertex->get_y();
0701     const auto dz = z - vertex->get_z();
0702     const auto pathlength = std::sqrt(square(dx) + square(dy) + square(dz));
0703     if (std::isfinite(pathlength) && pathlength > 0)
0704     {
0705       return pathlength;
0706     }
0707   }
0708 
0709   return static_cast<float>(stateIndex);
0710 }
0711 
0712 int PHTruthTrackFitter::getCharge(const PHG4Particle* particle,
0713                                   const TrackSeed* tpcSeed,
0714                                   const TrackSeed* siliconSeed) const
0715 {
0716   for (const auto* seed : {tpcSeed, siliconSeed})
0717   {
0718     if (!seed)
0719     {
0720       continue;
0721     }
0722 
0723     const auto charge = seed->get_charge();
0724     if (std::abs(charge) == 1)
0725     {
0726       return charge;
0727     }
0728   }
0729 
0730   const auto* pdg_particle = particle ? TDatabasePDG::Instance()->GetParticle(particle->get_pid()) : nullptr;
0731   if (pdg_particle && pdg_particle->Charge() < 0)
0732   {
0733     return -1;
0734   }
0735 
0736   return 1;
0737 }
0738 
0739 short int PHTruthTrackFitter::getCrossing(const TrackSeed* tpcSeed, const TrackSeed* siliconSeed) const
0740 {
0741   for (const auto* seed : {siliconSeed, tpcSeed})
0742   {
0743     if (!seed)
0744     {
0745       continue;
0746     }
0747 
0748     const auto crossing = seed->get_crossing();
0749     if (crossing != std::numeric_limits<short int>::max())
0750     {
0751       return crossing;
0752     }
0753   }
0754 
0755   return m_defaultCrossing;
0756 }
[ 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!