packages/trackreco/PHSiliconTrackMatching.cc

0001 #include "PHSiliconTrackMatching.h"
0002
0003 /// Tracking includes
0004 #include <trackbase/MvtxDefs.h>
0005 #include <trackbase/TrackFitUtils.h>
0006 #include <trackbase/TpcDefs.h>
0007 #include <trackbase/TrkrClusterContainer.h>
0008 #include <trackbase/TrkrClusterCrossingAssoc.h>
0009 #include <trackbase/TrkrClusterv3.h>
0010 #include <trackbase/TrkrDefs.h>  // for cluskey, getTrkrId, tpcId
0011
0012 #include <trackbase_historic/SvtxTrackSeed_v2.h>
0013 #include <trackbase_historic/TrackSeedContainer_v1.h>
0014 #include <trackbase_historic/TrackSeed_v2.h>
0015 #include <trackbase_historic/TrackSeedHelper.h>
0016
0017 #include <globalvertex/SvtxVertex.h>  // for SvtxVertex
0018 #include <globalvertex/SvtxVertexMap.h>
0019
0020 #include <g4main/PHG4Hit.h>       // for PHG4Hit
0021 #include <g4main/PHG4HitDefs.h>   // for keytype
0022 #include <g4main/PHG4Particle.h>  // for PHG4Particle
0023
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025
0026 #include <phool/PHCompositeNode.h>
0027 #include <phool/getClass.h>
0028 #include <phool/phool.h>
0029 #include <phool/sphenix_constants.h>
0030
0031 #include <TF1.h>
0032 #include <TFile.h>
0033 #include <TNtuple.h>
0034
0035 #include <climits>   // for UINT_MAX
0036 #include <cmath>     // for fabs, sqrt
0037 #include <iostream>  // for operator<<, basic_ostream
0038 #include <memory>
0039 #include <set>      // for _Rb_tree_const_iterator
0040 #include <utility>  // for pair
0041
0042 //____________________________________________________________________________..
0043 PHSiliconTrackMatching::PHSiliconTrackMatching(const std::string &name)
0044   : SubsysReco(name)
0045 {
0046 }
0047
0048 //____________________________________________________________________________..
0049 PHSiliconTrackMatching::~PHSiliconTrackMatching() = default;
0050
0051 //____________________________________________________________________________..
0052 int PHSiliconTrackMatching::InitRun(PHCompositeNode *topNode)
0053 {
0054   int ret = GetNodes(topNode);
0055   if (ret != Fun4AllReturnCodes::EVENT_OK)
0056   {
0057     return ret;
0058   }
0059   return ret;
0060 }
0061
0062
0063 //____________________________________________________________________________..
0064 int PHSiliconTrackMatching::process_event(PHCompositeNode * /*unused*/)
0065 {
0066   if(Verbosity() > 2)
0067   {
0068     std::cout << " Warning: PHSiliconTpcTrackMatching "
0069       << ( _zero_field ? "zero field is ON" : " zero field is OFF") << std::endl;
0070   }
0071   // _track_map contains the TPC seed track stubs
0072   // _track_map_silicon contains the silicon seed track stubs
0073   // _svtx_seed_map contains the combined silicon and tpc track seeds
0074
0075     // in case these objects are in the input file, we clear the nodes and replace them
0076   _svtx_seed_map->Reset();
0077   _track_map->Reset();
0078
0079   if (Verbosity() > 0)
0080   {
0081     std::cout << PHWHERE << " TPC track map size " << _track_map->size() << " Silicon track map size " << _track_map_silicon->size() << std::endl;
0082   }
0083
0084   if (_track_map_silicon->empty())
0085   {
0086     return Fun4AllReturnCodes::EVENT_OK;
0087   }
0088
0089   unsigned int si_id = 0;
0090   // loop over the silicon seeds and add the crossing to them
0091   for (unsigned int trackid = 0; trackid != _track_map_silicon->size(); ++trackid)
0092   {
0093     _tracklet_si = _track_map_silicon->get(trackid);
0094     if (!_tracklet_si)
0095     {
0096       continue;
0097     }
0098     auto crossing = _tracklet_si->get_crossing();
0099     if (Verbosity() > 8)
0100     {
0101       std::cout << " silicon stub: " << trackid << " eta " << _tracklet_si->get_eta()
0102         << " pt " << _tracklet_si->get_pt() << " si z " << TrackSeedHelper::get_z(_tracklet_si)
0103         << " crossing " << crossing << std::endl;
0104     }
0105
0106     if (Verbosity() > 1)
0107     {
0108       std::cout << " Si track " << trackid << " crossing " << crossing << std::endl;
0109     }
0110     auto dummy = std::make_unique<TrackSeed_v2>();
0111     dummy->set_qOverR(_tracklet_si->get_qOverR());
0112     dummy->set_phi(_tracklet_si->get_phi());
0113     dummy->set_X0(_tracklet_si->get_X0());
0114     dummy->set_Y0(_tracklet_si->get_Y0());
0115     dummy->set_Z0(_tracklet_si->get_Z0());
0116     dummy->set_slope(_tracklet_si->get_slope());
0117
0118     auto svtxseed = std::make_unique<SvtxTrackSeed_v2>();
0119     svtxseed->set_silicon_seed_index(si_id);
0120     svtxseed->set_tpc_seed_index(si_id);
0121     // In pp mode, if a matched track does not have INTT clusters we have to find the crossing geometrically
0122     // Record the geometrically estimated crossing in the track seeds for later use if needed
0123     svtxseed->set_crossing_estimate(crossing);
0124     _track_map->insert(dummy.get());
0125     _svtx_seed_map->insert(svtxseed.get());
0126
0127     si_id++;
0128     if (Verbosity() > 1)
0129     {
0130       std::cout << "  combined seed id " << _svtx_seed_map->size() - 1 << " si id " << si_id << " tpc id " << si_id << " crossing estimate " << crossing << std::endl;
0131     }
0132   }
0133
0134
0135   if (Verbosity() > 0)
0136   {
0137     std::cout << "final svtx seed map size " << _svtx_seed_map->size() << std::endl;
0138   }
0139
0140   if (Verbosity() > 1)
0141   {
0142     for (const auto &seed : *_svtx_seed_map)
0143     {
0144       seed->identify();
0145       std::cout << std::endl;
0146     }
0147
0148     std::cout << "PHSiliconTrackMatching::process_event(PHCompositeNode *topNode) Leaving process_event" << std::endl;
0149   }
0150   m_event++;
0151   return Fun4AllReturnCodes::EVENT_OK;
0152 }
0153
0154
0155 int PHSiliconTrackMatching::End(PHCompositeNode * /*unused*/)
0156 {
0157   return Fun4AllReturnCodes::EVENT_OK;
0158 }
0159
0160 int PHSiliconTrackMatching::GetNodes(PHCompositeNode *topNode)
0161 {
0162   //---------------------------------
0163   // Get additional objects off the Node Tree
0164   //---------------------------------
0165
0166   _track_map_silicon = findNode::getClass<TrackSeedContainer>(topNode, _silicon_track_map_name);
0167   if (!_track_map_silicon)
0168   {
0169     std::cerr << PHWHERE << " ERROR: Can't find SiliconTrackSeedContainer " << std::endl;
0170     return Fun4AllReturnCodes::ABORTEVENT;
0171   }
0172
0173   _track_map = findNode::getClass<TrackSeedContainer>(topNode, _track_map_name);
0174   if (!_track_map)
0175   {
0176     std::cerr << PHWHERE << " ERROR: Can't find " << _track_map_name << std::endl;
0177     return Fun4AllReturnCodes::ABORTEVENT;
0178   }
0179
0180   _svtx_seed_map = findNode::getClass<TrackSeedContainer>(topNode, "SvtxTrackSeedContainer");
0181   if (!_svtx_seed_map)
0182   {
0183     std::cout << "Creating node SvtxTrackSeedContainer" << std::endl;
0184     /// Get the DST Node
0185     PHNodeIterator iter(topNode);
0186     PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0187
0188     /// Check that it is there
0189     if (!dstNode)
0190     {
0191       std::cerr << "DST Node missing, quitting" << std::endl;
0192       throw std::runtime_error("failed to find DST node in PHActsSourceLinks::createNodes");
0193     }
0194
0195     /// Get the tracking subnode
0196     PHNodeIterator dstIter(dstNode);
0197     PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0198
0199     /// Check that it is there
0200     if (!svtxNode)
0201     {
0202       svtxNode = new PHCompositeNode("SVTX");
0203       dstNode->addNode(svtxNode);
0204     }
0205
0206     _svtx_seed_map = new TrackSeedContainer_v1();
0207     PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(_svtx_seed_map, "SvtxTrackSeedContainer", "PHObject");
0208     svtxNode->addNode(node);
0209   }
0210
0211   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, _cluster_map_name);
0212   if (!_cluster_map)
0213   {
0214     std::cout << PHWHERE << " ERROR: Can't find node " <<_cluster_map_name << std::endl;
0215     return Fun4AllReturnCodes::ABORTEVENT;
0216   }
0217
0218   _tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0219   if (!_tGeometry)
0220   {
0221     std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
0222     return Fun4AllReturnCodes::ABORTEVENT;
0223   }
0224
0225   return Fun4AllReturnCodes::EVENT_OK;
0226 }
0227