packages/trackreco/PHCosmicSiliconPropagator.cc

0001
0002 #include "PHCosmicSiliconPropagator.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/getClass.h>
0007 #include <trackbase/ActsGeometry.h>
0008 #include <trackbase/TrackFitUtils.h>
0009 #include <trackbase/TrkrCluster.h>
0010 #include <trackbase/TrkrClusterContainer.h>
0011 #include <trackbase/TrkrClusterCrossingAssoc.h>
0012 #include <trackbase_historic/SvtxTrackMap.h>
0013 #include <trackbase_historic/SvtxTrackSeed_v1.h>
0014 #include <trackbase_historic/TrackSeedContainer.h>
0015 #include <trackbase_historic/TrackSeedContainer_v1.h>
0016 #include <trackbase_historic/TrackSeed_v2.h>
0017 #include <trackbase_historic/TrackSeedHelper.h>
0018
0019 #include <cmath>
0020
0021 namespace
0022 {
0023   template <class T>
0024   inline constexpr T square(T& x)
0025   {
0026     return x * x;
0027   }
0028   template <class T>
0029   inline constexpr T r(T& x, T& y)
0030   {
0031     return std::sqrt(square(x) + square(y));
0032   }
0033
0034 }  // namespace
0035 //____________________________________________________________________________..
0036 PHCosmicSiliconPropagator::PHCosmicSiliconPropagator(const std::string& name)
0037   : SubsysReco(name)
0038 {
0039 }
0040
0041 //____________________________________________________________________________..
0042 PHCosmicSiliconPropagator::~PHCosmicSiliconPropagator() = default;
0043
0044 //____________________________________________________________________________..
0045 int PHCosmicSiliconPropagator::Init(PHCompositeNode* /*unused*/)
0046 {
0047   return Fun4AllReturnCodes::EVENT_OK;
0048 }
0049
0050 //____________________________________________________________________________..
0051 int PHCosmicSiliconPropagator::InitRun(PHCompositeNode* topNode)
0052 {
0053   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0054   if (!_cluster_map)
0055   {
0056     std::cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
0057     return Fun4AllReturnCodes::ABORTEVENT;
0058   }
0059
0060   _tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0061   if (!_tgeometry)
0062   {
0063     std::cout << "No Acts tracking geometry, exiting." << std::endl;
0064     return Fun4AllReturnCodes::ABORTEVENT;
0065   }
0066   _tpc_seeds = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0067   if (!_tpc_seeds)
0068   {
0069     std::cout << "No TpcTrackSeedContainer, exiting." << std::endl;
0070     return Fun4AllReturnCodes::ABORTEVENT;
0071   }
0072   _si_seeds = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0073   if (!_si_seeds)
0074   {
0075     std::cout << "No SiliconTrackSeedContainer, creating..." << std::endl;
0076     if (createSeedContainer(_si_seeds, "SiliconTrackSeedContainer", topNode) != Fun4AllReturnCodes::EVENT_OK)
0077     {
0078       std::cout << "Cannot create, exiting." << std::endl;
0079       return Fun4AllReturnCodes::ABORTEVENT;
0080     }
0081   }
0082   _svtx_seeds = findNode::getClass<TrackSeedContainer>(topNode, _track_map_name);
0083   if (!_svtx_seeds)
0084   {
0085     std::cout << "No " << _track_map_name << " found, creating..." << std::endl;
0086     if (createSeedContainer(_svtx_seeds, _track_map_name, topNode) != Fun4AllReturnCodes::EVENT_OK)
0087     {
0088       std::cout << "Cannot create, exiting." << std::endl;
0089       return Fun4AllReturnCodes::ABORTEVENT;
0090     }
0091   }
0092
0093   return Fun4AllReturnCodes::EVENT_OK;
0094 }
0095
0096 //____________________________________________________________________________..
0097 int PHCosmicSiliconPropagator::process_event(PHCompositeNode* /*unused*/)
0098 {
0099   if (m_resetContainer)
0100   {
0101     _svtx_seeds->Reset();
0102   }
0103
0104   for (auto& tpcseed : *_tpc_seeds)
0105   {
0106     if (!tpcseed)
0107     {
0108       continue;
0109     }
0110     if (Verbosity() > 3)
0111     {
0112       std::cout << "Tpc seed to start out is " << std::endl;
0113       tpcseed->identify();
0114     }
0115     std::vector<Acts::Vector3> tpcClusPos;
0116     std::vector<TrkrDefs::cluskey> tpcClusKeys;
0117
0118     std::copy(tpcseed->begin_cluster_keys(), tpcseed->end_cluster_keys(),
0119               std::back_inserter(tpcClusKeys));
0120
0121     TrackFitUtils::getTrackletClusters(_tgeometry, _cluster_map,
0122                                        tpcClusPos, tpcClusKeys);
0123     std::vector<TrkrDefs::cluskey> newClusKeys;
0124     std::vector<Acts::Vector3> newClusPos;
0125     unsigned int nClusters = 0;
0126
0127     if (m_zeroField)
0128     {
0129       //! Line fit both in x-y and r-z
0130       TrackFitUtils::position_vector_t rzpoints, xypoints;
0131       for (auto& globPos : tpcClusPos)
0132       {
0133         xypoints.push_back(std::make_pair(globPos.x(), globPos.y()));
0134         float clusr = r(globPos.x(), globPos.y());
0135         if (globPos.y() < 0)
0136         {
0137           clusr *= -1;
0138         }
0139         rzpoints.push_back(std::make_pair(globPos.z(), clusr));
0140       }
0141
0142       auto xyLineParams = TrackFitUtils::line_fit(xypoints);
0143       auto rzLineParams = TrackFitUtils::line_fit(rzpoints);
0144
0145       std::vector<TrkrDefs::cluskey> newClusKeysrz;
0146       std::vector<Acts::Vector3> newClusPosrz;
0147       std::vector<TrkrDefs::cluskey> newClusKeysxy;
0148       std::vector<Acts::Vector3> newClusPosxy;
0149       // now add clusters along lines
0150       nClusters = TrackFitUtils::addClustersOnLine(xyLineParams,
0151                                                    true,
0152                                                    _dca_xy_cut,
0153                                                    _tgeometry,
0154                                                    _cluster_map,
0155                                                    newClusPosxy,
0156                                                    newClusKeysxy,
0157                                                    0, 56);
0158       int nrzClusters = TrackFitUtils::addClustersOnLine(rzLineParams,
0159                                                          false,
0160                                                          _dca_z_cut,
0161                                                          _tgeometry,
0162                                                          _cluster_map,
0163                                                          newClusPosrz,
0164                                                          newClusKeysrz,
0165                                                          0, 56);
0166
0167       if (Verbosity() > 3)
0168       {
0169         std::cout << "nrz clus " << nrzClusters << " and nxy clusters " << nClusters
0170                   << std::endl;
0171       }
0172       std::set_intersection(newClusKeysxy.begin(), newClusKeysxy.end(),
0173                             newClusKeysrz.begin(), newClusKeysrz.end(), std::back_inserter(newClusKeys));
0174       if (m_resetContainer)
0175       {
0176         for (auto& keys : {newClusKeysxy, newClusKeysrz})
0177         {
0178           for (auto& key : keys)
0179           {
0180             if (TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::micromegasId)
0181             {
0182               newClusKeys.push_back(key);
0183             }
0184           }
0185         }
0186       }
0187
0188       if (Verbosity() > 3)
0189       {
0190         for (auto key : newClusKeysxy)
0191         {
0192           auto cluster = _cluster_map->findCluster(key);
0193           auto clusglob = _tgeometry->getGlobalPosition(key, cluster);
0194           std::cout << "Found key " << key << " for xy cosmic in layer " << (unsigned int) TrkrDefs::getLayer(key)
0195                     << " with pos " << clusglob.transpose() << std::endl;
0196         }
0197         for (auto key : newClusKeysrz)
0198         {
0199           auto cluster = _cluster_map->findCluster(key);
0200           auto clusglob = _tgeometry->getGlobalPosition(key, cluster);
0201           std::cout << "Found key " << key << " for rz cosmic in layer " << (unsigned int) TrkrDefs::getLayer(key)
0202                     << " with pos " << clusglob.transpose() << std::endl;
0203         }
0204       }
0205     }
0206     else
0207     {
0208       std::vector<float> fitparams = TrackFitUtils::fitClusters(tpcClusPos, tpcClusKeys);
0209       //! There weren't enough clusters to fit
0210       if (fitparams.size() == 0)
0211       {
0212         continue;
0213       }
0214
0215       std::vector<TrkrDefs::cluskey> ckeys;
0216       nClusters = TrackFitUtils::addClusters(fitparams, _dca_xy_cut, _tgeometry, _cluster_map,
0217                                              newClusPos, ckeys, 0, 56);
0218       TrackFitUtils::position_vector_t yzpoints;
0219       for (auto& globPos : tpcClusPos)
0220       {
0221         yzpoints.push_back(std::make_pair(globPos.y(), globPos.z()));
0222       }
0223
0224       auto yzLineParams = TrackFitUtils::line_fit(yzpoints);
0225       float yzslope = std::get<0>(yzLineParams);
0226       float yzint = std::get<1>(yzLineParams);
0227       for (auto& key : ckeys)
0228       {
0229         auto cluster = _cluster_map->findCluster(key);
0230         auto clusglob = _tgeometry->getGlobalPosition(key, cluster);
0231
0232         float projz = clusglob.y() * yzslope + yzint;
0233
0234         if (std::fabs(projz - clusglob.z()) < _dca_z_cut)
0235         {
0236           newClusKeys.push_back(key);
0237         }
0238       }
0239     }
0240     //! only keep long seeds
0241     if ((tpcClusKeys.size() + newClusKeys.size() > 25))
0242     {
0243       // TODO: should include distortion corrections
0244       std::unique_ptr<TrackSeed_v2> si_seed = std::make_unique<TrackSeed_v2>();
0245       std::map<TrkrDefs::cluskey, Acts::Vector3> silposmap, tpcposmap;
0246       for (auto& key : tpcClusKeys)
0247       {
0248         auto cluster = _cluster_map->findCluster(key);
0249         auto clusglob = _tgeometry->getGlobalPosition(key, cluster);
0250         tpcposmap.emplace(key, clusglob);
0251       }
0252       for (auto& key : newClusKeys)
0253       {
0254         bool isTpcKey = false;
0255         auto cluster = _cluster_map->findCluster(key);
0256         auto clusglob = _tgeometry->getGlobalPosition(key, cluster);
0257         if (TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::tpcId ||
0258             TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::micromegasId)
0259         {
0260           isTpcKey = true;
0261         }
0262         if (!isTpcKey)
0263         {
0264           si_seed->insert_cluster_key(key);
0265           silposmap.emplace(key, clusglob);
0266         }
0267         else
0268         {
0269           tpcseed->insert_cluster_key(key);
0270           tpcposmap.emplace(key, clusglob);
0271         }
0272       }
0273
0274       TrackSeedHelper::circleFitByTaubin(si_seed.get(), silposmap, 0, 8);
0275       TrackSeedHelper::lineFit(si_seed.get(), silposmap);
0276
0277       TrackSeedHelper::circleFitByTaubin(tpcseed, tpcposmap, 7, 57);
0278       TrackSeedHelper::lineFit(tpcseed, tpcposmap, 7, 57);
0279
0280       TrackSeed* mapped_seed = _si_seeds->insert(si_seed.get());
0281       std::unique_ptr<SvtxTrackSeed_v1> full_seed = std::make_unique<SvtxTrackSeed_v1>();
0282       int tpcind = _tpc_seeds->find(tpcseed);
0283       int siind = _si_seeds->find(mapped_seed);
0284       full_seed->set_tpc_seed_index(tpcind);
0285       if (si_seed->size_cluster_keys() > 0)
0286       {
0287         full_seed->set_silicon_seed_index(siind);
0288       }
0289       _svtx_seeds->insert(full_seed.get());
0290       if (Verbosity() > 3)
0291       {
0292         std::cout << "final seeds" << std::endl;
0293         si_seed->identify();
0294         tpcseed->identify();
0295       }
0296     }
0297   }
0298
0299   if (Verbosity() > 2)
0300   {
0301     std::cout << "svtx seed map size is " << _svtx_seeds->size() << std::endl;
0302     int i = 0;
0303     for (auto& seed : *_svtx_seeds)
0304     {
0305       std::cout << "seed " << i << " is composed of " << std::endl;
0306       _tpc_seeds->get(seed->get_tpc_seed_index())->identify();
0307       if (_si_seeds->get(seed->get_silicon_seed_index()))
0308       {
0309         _si_seeds->get(seed->get_silicon_seed_index())->identify();
0310       }
0311       ++i;
0312     }
0313   }
0314   return Fun4AllReturnCodes::EVENT_OK;
0315 }
0316
0317 //____________________________________________________________________________..
0318 int PHCosmicSiliconPropagator::End(PHCompositeNode* /*unused*/)
0319 {
0320   return Fun4AllReturnCodes::EVENT_OK;
0321 }
0322
0323 int PHCosmicSiliconPropagator::createSeedContainer(TrackSeedContainer*& container, const std::string& container_name, PHCompositeNode* topNode)
0324 {
0325   PHNodeIterator iter(topNode);
0326
0327   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0328
0329   if (!dstNode)
0330   {
0331     std::cerr << "DST node is missing, quitting" << std::endl;
0332     throw std::runtime_error("Failed to find DST node in PHSiliconHelicalPropagator::createNodes");
0333   }
0334
0335   PHNodeIterator dstIter(dstNode);
0336   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0337
0338   if (!svtxNode)
0339   {
0340     svtxNode = new PHCompositeNode("SVTX");
0341     dstNode->addNode(svtxNode);
0342   }
0343
0344   container = findNode::getClass<TrackSeedContainer>(topNode, container_name);
0345   if (!container)
0346   {
0347     container = new TrackSeedContainer_v1;
0348     PHIODataNode<PHObject>* trackNode = new PHIODataNode<PHObject>(container, container_name, "PHObject");
0349     svtxNode->addNode(trackNode);
0350   }
0351
0352   return Fun4AllReturnCodes::EVENT_OK;
0353 }