packages/trackreco/PHSiliconHelicalPropagator.cc

0001 #include "PHSiliconHelicalPropagator.h"
0002
0003 #include <trackbase/TrkrClusterContainer.h>
0004 #include <trackbase_historic/SvtxTrackSeed_v1.h>
0005 #include <trackbase_historic/TrackSeedContainer_v1.h>
0006 #include <trackbase_historic/TrackSeed_v2.h>
0007 #include <trackbase_historic/TrackSeedHelper.h>
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <phool/PHCompositeNode.h>
0011 #include <phool/getClass.h>
0012
0013 namespace
0014 {
0015   template <typename T>
0016   int sgn(const T& x)
0017   {
0018     if (x > 0)
0019     {
0020       return 1;
0021     }
0022     else
0023     {
0024       return -1;
0025     }
0026   }
0027 }  // namespace
0028
0029 PHSiliconHelicalPropagator::PHSiliconHelicalPropagator(const std::string& name)
0030   : SubsysReco(name)
0031 {
0032 }
0033
0034 PHSiliconHelicalPropagator::~PHSiliconHelicalPropagator() = default;
0035
0036 int PHSiliconHelicalPropagator::InitRun(PHCompositeNode* topNode)
0037 {
0038   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0039   if (!_cluster_map)
0040   {
0041     std::cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
0042     return Fun4AllReturnCodes::ABORTEVENT;
0043   }
0044   _cluster_crossing_map = findNode::getClass<TrkrClusterCrossingAssoc>(topNode, "TRKR_CLUSTERCROSSINGASSOC");
0045   if (!_cluster_crossing_map)
0046   {
0047     std::cerr << PHWHERE << " ERROR: Can't find TRKR_CLUSTERCROSSINGASSOC " << std::endl;
0048     return Fun4AllReturnCodes::ABORTEVENT;
0049   }
0050   _tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0051   if (!_tgeometry)
0052   {
0053     std::cout << "No Acts tracking geometry, exiting." << std::endl;
0054     return Fun4AllReturnCodes::ABORTEVENT;
0055   }
0056   _tpc_seeds = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0057   if (!_tpc_seeds)
0058   {
0059     std::cout << "No TpcTrackSeedContainer, exiting." << std::endl;
0060     return Fun4AllReturnCodes::ABORTEVENT;
0061   }
0062   _si_seeds = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0063   if (!_si_seeds)
0064   {
0065     std::cout << "No SiliconTrackSeedContainer, creating..." << std::endl;
0066     if (createSeedContainer(_si_seeds, "SiliconTrackSeedContainer", topNode) != Fun4AllReturnCodes::EVENT_OK)
0067     {
0068       std::cout << "Cannot create, exiting." << std::endl;
0069       return Fun4AllReturnCodes::ABORTEVENT;
0070     }
0071   }
0072   _svtx_seeds = findNode::getClass<TrackSeedContainer>(topNode, _track_map_name);
0073   if (!_svtx_seeds)
0074   {
0075     std::cout << "No " << _track_map_name << " found, creating..." << std::endl;
0076     if (createSeedContainer(_svtx_seeds, _track_map_name, topNode) != Fun4AllReturnCodes::EVENT_OK)
0077     {
0078       std::cout << "Cannot create, exiting." << std::endl;
0079       return Fun4AllReturnCodes::ABORTEVENT;
0080     }
0081   }
0082   return Fun4AllReturnCodes::EVENT_OK;
0083 }
0084
0085 int PHSiliconHelicalPropagator::createSeedContainer(TrackSeedContainer*& container, const std::string& container_name, PHCompositeNode* topNode)
0086 {
0087   PHNodeIterator iter(topNode);
0088
0089   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0090
0091   if (!dstNode)
0092   {
0093     std::cerr << "DST node is missing, quitting" << std::endl;
0094     throw std::runtime_error("Failed to find DST node in PHSiliconHelicalPropagator::createNodes");
0095   }
0096
0097   PHNodeIterator dstIter(dstNode);
0098   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0099
0100   if (!svtxNode)
0101   {
0102     svtxNode = new PHCompositeNode("SVTX");
0103     dstNode->addNode(svtxNode);
0104   }
0105
0106   container = findNode::getClass<TrackSeedContainer>(topNode, container_name);
0107   if (!container)
0108   {
0109     container = new TrackSeedContainer_v1;
0110     PHIODataNode<PHObject>* trackNode = new PHIODataNode<PHObject>(container, container_name, "PHObject");
0111     svtxNode->addNode(trackNode);
0112   }
0113
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116
0117 int PHSiliconHelicalPropagator::process_event(PHCompositeNode* /*topNode*/)
0118 {
0119   for (unsigned int seedID = 0; seedID < _tpc_seeds->size(); ++seedID)
0120   {
0121     TrackSeed* tpcseed = _tpc_seeds->get(seedID);
0122     if (!tpcseed)
0123     {
0124       continue;
0125     }
0126   if(Verbosity() > 2)
0127   {
0128     tpcseed->identify();
0129   }
0130     std::vector<Acts::Vector3> clusterPositions;
0131     std::vector<TrkrDefs::cluskey> clusterKeys;
0132     for (auto iter = tpcseed->begin_cluster_keys();
0133          iter != tpcseed->end_cluster_keys(); ++iter)
0134     {
0135       clusterKeys.push_back(*iter);
0136     }
0137     TrackFitUtils::getTrackletClusters(_tgeometry, _cluster_map, clusterPositions, clusterKeys);
0138     std::vector<float> fitparams = TrackFitUtils::fitClusters(clusterPositions, clusterKeys);
0139     //! There weren't enough clusters to fit
0140     if (fitparams.size() == 0)
0141     {
0142       continue;
0143     }
0144     if(Verbosity() > 3)
0145     {
0146       for(auto& param : fitparams)
0147       {
0148         std::cout << "fit param " << param << std::endl;
0149       }
0150     }
0151     std::vector<TrkrDefs::cluskey> si_clusterKeys, si_clusterKeysrz;
0152     std::vector<Acts::Vector3> si_clusterPositions, si_clusterPositionsrz;
0153     std::map<TrkrDefs::cluskey, Acts::Vector3> positionMap;
0154
0155     unsigned int nSiClusters = std::numeric_limits<unsigned int>::quiet_NaN();
0156     TrackFitUtils::position_vector_t xypoints, rzpoints;
0157     for (auto& pos : clusterPositions)
0158     {
0159       xypoints.push_back({pos.x(), pos.y()});
0160       float clusr = std::sqrt(pos.x() * pos.x() + pos.y() * pos.y());
0161       if (pos.y() < 0)
0162       {
0163         clusr = -clusr;
0164       }
0165       rzpoints.push_back({pos.z(), clusr});
0166
0167       if (Verbosity() > 5)
0168       {
0169         std::cout << "Cluster pos " << pos.transpose() << " and r " << std::sqrt(pos.x() * pos.x() + pos.y() * pos.y()) << std::endl;
0170       }
0171     }
0172     auto rzparams = TrackFitUtils::line_fit(rzpoints);
0173
0174     if (m_zeroField)
0175     {
0176
0177      auto xyparams = TrackFitUtils::line_fit(xypoints);
0178      nSiClusters = TrackFitUtils::addClustersOnLine(xyparams,
0179                                                     true,
0180                                                     _dca_cut,
0181                                                     _tgeometry, _cluster_map,
0182                                                     si_clusterPositions,
0183                                                     si_clusterKeys,
0184                                                     0, 6);
0185    }
0186    else{
0187     nSiClusters = TrackFitUtils::addClusters(fitparams, _dca_cut, _tgeometry, _cluster_map, si_clusterPositions, si_clusterKeys, 0, 6);
0188    }
0189    int nrzClusters = TrackFitUtils::addClustersOnLine(rzparams,
0190                                                        false,
0191                                                        _dca_z_cut,
0192                                                        _tgeometry,
0193                                                        _cluster_map,
0194                                                        si_clusterPositionsrz,
0195                                                        si_clusterKeysrz,
0196                                                        0, 6);
0197     std::vector<TrkrDefs::cluskey> newkeys;
0198     std::set_intersection(si_clusterKeys.begin(), si_clusterKeys.end(),
0199                           si_clusterKeysrz.begin(), si_clusterKeysrz.end(),
0200                           std::back_inserter(newkeys));
0201
0202     if (newkeys.size() > 0)
0203     {
0204       if(Verbosity() > 0)
0205       {
0206         std::cout << "Adding " << newkeys.size() << " Keys " << std::endl;
0207         for(auto& key : newkeys)
0208         {
0209           std::cout << "key " << (unsigned int) key << std::endl;
0210         }
0211       }
0212       std::unique_ptr<TrackSeed_v2> si_seed = std::make_unique<TrackSeed_v2>();
0213       std::map<short, int> crossing_frequency;
0214       Acts::Vector3 tpcExGlobal = clusterPositions.front();
0215       for (auto& clusterkey : newkeys)
0216       {
0217         //! Check that the clusters are in the same quadrant
0218         auto cluster = _cluster_map->findCluster(clusterkey);
0219         auto global = _tgeometry->getGlobalPosition(clusterkey, cluster);
0220         positionMap.insert({clusterkey, global});
0221         if (sgn(global.x()) == sgn(tpcExGlobal.x()) && sgn(global.y()) == sgn(tpcExGlobal.y()))
0222         {
0223           si_seed->insert_cluster_key(clusterkey);
0224         }
0225         /*
0226         else if (TrkrDefs::getTrkrId(clusterkey) == TrkrDefs::inttId)
0227         {
0228           auto hit_crossings = _cluster_crossing_map->getCrossings(clusterkey);
0229           for (auto iter = hit_crossings.first; iter != hit_crossings.second; ++iter)
0230           {
0231             short crossing = iter->second;
0232             if (crossing_frequency.count(crossing) == 0)
0233             {
0234               crossing_frequency.insert({crossing, 1});
0235             }
0236             else
0237             {
0238               crossing_frequency[crossing]++;
0239             }
0240           }
0241
0242
0243           if (sgn(global.x()) == sgn(tpcExGlobal.x()) && sgn(global.y()) == sgn(tpcExGlobal.y()))
0244           {
0245             si_seed->insert_cluster_key(clusterkey);
0246           }
0247         }
0248         */
0249       }
0250 /*
0251       if (crossing_frequency.size() > 0)
0252       {
0253         short most_common_crossing = (std::max_element(crossing_frequency.begin(), crossing_frequency.end(),
0254                                                        [](auto entry1, auto entry2)
0255                                                        { return entry1.second > entry2.second; }))
0256                                          ->first;
0257         si_seed->set_crossing(most_common_crossing);
0258       }
0259       */
0260       TrackSeedHelper::circleFitByTaubin(si_seed.get(), positionMap, 0, 8);
0261       TrackSeedHelper::lineFit(si_seed.get(), positionMap, 0, 8);
0262       si_seed->set_crossing(0);
0263       TrackSeed* mapped_seed = _si_seeds->insert(si_seed.get());
0264
0265       std::unique_ptr<SvtxTrackSeed_v1> full_seed = std::make_unique<SvtxTrackSeed_v1>();
0266       int tpc_seed_index = _tpc_seeds->find(tpcseed);
0267       int si_seed_index = _si_seeds->find(mapped_seed);
0268       if (Verbosity() > 0)
0269       {
0270         std::cout << "found  " << nSiClusters << " silicon clusters in xy for tpc seed " << tpc_seed_index << std::endl;
0271         std::cout << "found " << nrzClusters << " silicon clusters in rz for tpc seed " << tpc_seed_index << std::endl;
0272         std::cout << "intersection is " << newkeys.size() << std::endl;
0273         std::cout << "new silicon seed index: " << si_seed_index << std::endl;
0274       }
0275       full_seed->set_tpc_seed_index(tpc_seed_index);
0276       full_seed->set_silicon_seed_index(si_seed_index);
0277       _svtx_seeds->insert(full_seed.get());
0278     }
0279     else
0280     {
0281       // no Si clusters found, put TPC-only seed in SvtxTrackSeedContainer
0282       std::unique_ptr<SvtxTrackSeed_v1> partial_seed = std::make_unique<SvtxTrackSeed_v1>();
0283       int tpc_seed_index = _tpc_seeds->find(tpcseed);
0284       partial_seed->set_tpc_seed_index(tpc_seed_index);
0285       _svtx_seeds->insert(partial_seed.get());
0286     }
0287   }
0288   if (Verbosity() > 2)
0289   {
0290     std::cout << "svtx seed map size is " << _svtx_seeds->size() << std::endl;
0291   }
0292   return Fun4AllReturnCodes::EVENT_OK;
0293 }
0294
0295 int PHSiliconHelicalPropagator::End(PHCompositeNode* /*topNode*/)
0296 {
0297   return Fun4AllReturnCodes::EVENT_OK;
0298 }