packages/trackreco/PHTrackPruner.cc

0001 #include "PHTrackPruner.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/TrkrClusterv3.h>
0009 #include <trackbase/TrkrDefs.h>  // for cluskey, getTrkrId, tpcId
0010
0011 #include <trackbase_historic/SvtxTrackSeed_v2.h>
0012 #include <trackbase_historic/TrackSeedContainer_v1.h>
0013 #include <trackbase_historic/TrackSeed_v2.h>
0014 #include <trackbase_historic/TrackSeedHelper.h>
0015
0016 #include <trackbase_historic/SvtxTrack.h>  // for SvtxTrack
0017 #include <trackbase_historic/SvtxTrackMap.h>
0018
0019 #include <fun4all/Fun4AllReturnCodes.h>
0020
0021 #include <phool/PHCompositeNode.h>
0022 #include <phool/getClass.h>
0023 #include <phool/phool.h>
0024 #include <phool/sphenix_constants.h>
0025
0026 #include <TF1.h>
0027 #include <TFile.h>
0028 #include <TNtuple.h>
0029
0030 #include <climits>   // for UINT_MAX
0031 #include <cmath>     // for fabs, sqrt
0032 #include <iostream>  // for operator<<, basic_ostream
0033 #include <memory>
0034 #include <set>      // for _Rb_tree_const_iterator
0035 #include <utility>  // for pair
0036
0037 using namespace std;
0038
0039 namespace
0040 {
0041   //! get cluster keys from a given track
0042   std::vector<TrkrDefs::cluskey> get_cluster_keys(SvtxTrack* track)
0043   {
0044     std::vector<TrkrDefs::cluskey> out;
0045     for (const auto& seed : {track->get_silicon_seed(), track->get_tpc_seed()})
0046     {
0047       if (seed)
0048       {
0049         std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
0050       }
0051     }
0052     return out;
0053   }
0054
0055   //! get cluster keys from a given track associated with states
0056   std::vector<TrkrDefs::cluskey> get_state_keys(SvtxTrack* track)
0057   {
0058     std::vector<TrkrDefs::cluskey> out;
0059     for (auto state_iter = track->begin_states();
0060          state_iter != track->end_states();
0061          ++state_iter)
0062     {
0063       SvtxTrackState* tstate = state_iter->second;
0064       auto stateckey = tstate->get_cluskey();
0065       out.push_back(stateckey);
0066     }
0067     return out;
0068   }
0069
0070   /// return number of clusters of a given type that belong to a tracks
0071   template <int type>
0072   int count_clusters(const std::vector<TrkrDefs::cluskey>& keys)
0073   {
0074     return std::count_if(keys.begin(), keys.end(),
0075                          [](const TrkrDefs::cluskey& key)
0076                          { return TrkrDefs::getTrkrId(key) == type; });
0077   }
0078 }
0079
0080 //____________________________________________________________________________..
0081 PHTrackPruner::PHTrackPruner(const std::string &name)
0082   : SubsysReco(name)
0083 {
0084 }
0085
0086 //____________________________________________________________________________..
0087 PHTrackPruner::~PHTrackPruner() = default;
0088
0089 //____________________________________________________________________________..
0090 int PHTrackPruner::InitRun(PHCompositeNode *topNode)
0091 {
0092   int ret = GetNodes(topNode);
0093   if (ret != Fun4AllReturnCodes::EVENT_OK)
0094   {
0095     return ret;
0096   }
0097
0098   return ret;
0099 }
0100
0101 //____________________________________________________________________________..
0102 int PHTrackPruner::process_event(PHCompositeNode * /*unused*/)
0103 {
0104   // _tpc_seed_map contains the TPC seed track stubs
0105   // _si_seed_map contains the silicon seed track stubs
0106   // _svtx_seed_map contains the combined silicon and tpc track seeds
0107   // _svtx_track_map contains the fitted acts track stubs
0108
0109   if (Verbosity() > 0)
0110   {
0111     cout << PHWHERE << " TPC seed map size " << _tpc_seed_map->size()
0112      << " Silicon seed map size " << _si_seed_map->size()
0113      << " Svtx seed map size " << _svtx_seed_map->size()
0114      << " Svtx track map size " << _svtx_track_map->size()
0115      << endl;
0116   }
0117
0118   if (_svtx_track_map->size() == 0)
0119   {
0120     return Fun4AllReturnCodes::EVENT_OK;
0121   }
0122
0123   std::multimap<unsigned int, unsigned int> good_matches;
0124
0125   for (auto &iter : *_svtx_track_map)
0126   {
0127     _svtx_track = iter.second;
0128
0129     if(!checkTrack(_svtx_track))
0130     {
0131       continue;
0132     }
0133     if (Verbosity() > 1) { std::cout<<"Pass track selection"<<std::endl; }
0134
0135     _tpc_seed = _svtx_track->get_tpc_seed();
0136     _si_seed = _svtx_track->get_silicon_seed();
0137     if (_tpc_seed && _si_seed)
0138     {
0139       if (Verbosity() > 1) { std::cout<<"Insert tpcid and siid into good_matches"<<std::endl; }
0140       int tpcid = _tpc_seed_map->find(_tpc_seed);
0141       int siid = _si_seed_map->find(_si_seed);
0142       good_matches.insert(std::make_pair(tpcid, siid));
0143     }
0144   }
0145
0146   for (auto [tpcid, siid] : good_matches)
0147   {
0148       if (Verbosity() > 1) { std::cout<<"Insert pruned svtx seed map"<<std::endl; }
0149     auto _svtx_seed = std::make_unique<SvtxTrackSeed_v2>();
0150     _svtx_seed->set_silicon_seed_index(siid);
0151     _svtx_seed->set_tpc_seed_index(tpcid);
0152     // In pp mode, if a matched track does not have INTT clusters we have to find the crossing geometrically
0153     // Record the geometrically estimated crossing in the track seeds for later use if needed
0154     short int crossing_estimate = findCrossingGeometrically(tpcid, siid);
0155     _svtx_seed->set_crossing_estimate(crossing_estimate);
0156     _pruned_svtx_seed_map->insert(_svtx_seed.get());
0157
0158     if (Verbosity() > 1)
0159     {
0160       std::cout << "  combined seed id " << _pruned_svtx_seed_map->size() - 1 << " si id " << siid << " tpc id " << tpcid << " crossing estimate " << crossing_estimate << std::endl;
0161     }
0162   }
0163
0164   if (Verbosity() > 0)
0165   {
0166     std::cout << "final svtx seed map size " << _pruned_svtx_seed_map->size() << std::endl;
0167   }
0168
0169   if (Verbosity() > 1)
0170   {
0171     for (const auto &seed : *_pruned_svtx_seed_map)
0172     {
0173       seed->identify();
0174     }
0175
0176     cout << "PHTrackPruner::process_event(PHCompositeNode *topNode) Leaving process_event" << endl;
0177   }
0178   m_event++;
0179   return Fun4AllReturnCodes::EVENT_OK;
0180 }
0181
0182 bool PHTrackPruner::checkTrack(SvtxTrack *track)
0183 {
0184   if(!track)
0185   {
0186     if (Verbosity() > 1) { std::cout<<"invalid track"<<std::endl; }
0187     return false;
0188   }
0189
0190   //pt cut
0191   if(track->get_pt() < m_track_pt_low_cut)
0192   {
0193     if (Verbosity() > 1) { std::cout<<"Track pt "<<track->get_pt()<<" , pt cut "<<m_track_pt_low_cut<<std::endl; }
0194     return false;
0195   }
0196
0197   //quality cut
0198   if(track->get_quality() > m_track_quality_high_cut)
0199   {
0200     if (Verbosity() > 1) { std::cout<<"Track quality "<<track->get_quality()<<" , quality cut "<<m_track_quality_high_cut<<std::endl; }
0201     return false;
0202   }
0203
0204   //number of clusters cut
0205   const auto cluster_keys(get_cluster_keys(track));
0206   if (count_clusters<TrkrDefs::mvtxId>(cluster_keys) < m_nmvtx_clus_low_cut)
0207   {
0208     if (Verbosity() > 1) { std::cout<<"nmvtx "<<count_clusters<TrkrDefs::mvtxId>(cluster_keys)<<" , nmvtx cut "<<m_nmvtx_clus_low_cut<<std::endl; }
0209     return false;
0210   }
0211   if (count_clusters<TrkrDefs::inttId>(cluster_keys) < m_nintt_clus_low_cut)
0212   {
0213     if (Verbosity() > 1) { std::cout<<"nintt "<<count_clusters<TrkrDefs::inttId>(cluster_keys)<<" , nintt cut "<<m_nintt_clus_low_cut<<std::endl; }
0214     return false;
0215   }
0216   if (count_clusters<TrkrDefs::tpcId>(cluster_keys) < m_ntpc_clus_low_cut)
0217   {
0218     if (Verbosity() > 1) { std::cout<<"ntpc "<<count_clusters<TrkrDefs::tpcId>(cluster_keys)<<" , ntpc cut "<<m_ntpc_clus_low_cut<<std::endl; }
0219     return false;
0220   }
0221   if (count_clusters<TrkrDefs::micromegasId>(cluster_keys) < m_ntpot_clus_low_cut)
0222   {
0223     if (Verbosity() > 1) { std::cout<<"nmicromegas "<<count_clusters<TrkrDefs::micromegasId>(cluster_keys)<<" , nmicromegas cut "<<m_ntpot_clus_low_cut<<std::endl; }
0224     return false;
0225   }
0226
0227   //number of states cut
0228   const auto state_keys(get_state_keys(track));
0229   if (count_clusters<TrkrDefs::mvtxId>(state_keys) < m_nmvtx_states_low_cut)
0230   {
0231     if (Verbosity() > 1) { std::cout<<"nmvtxstates "<<count_clusters<TrkrDefs::mvtxId>(state_keys)<<" , nmvtxstates cut "<<m_nmvtx_states_low_cut<<std::endl; }
0232     return false;
0233   }
0234   if (count_clusters<TrkrDefs::inttId>(state_keys) < m_nintt_states_low_cut)
0235   {
0236     if (Verbosity() > 1) { std::cout<<"ninttstates "<<count_clusters<TrkrDefs::inttId>(state_keys)<<" , ninttstates cut "<<m_nintt_states_low_cut<<std::endl; }
0237     return false;
0238   }
0239   if (count_clusters<TrkrDefs::tpcId>(state_keys) < m_ntpc_states_low_cut)
0240   {
0241     if (Verbosity() > 1) { std::cout<<"ntpcstates "<<count_clusters<TrkrDefs::tpcId>(state_keys)<<" , ntpcstates cut "<<m_ntpc_states_low_cut<<std::endl; }
0242     return false;
0243   }
0244   if (count_clusters<TrkrDefs::micromegasId>(state_keys) < m_ntpot_states_low_cut)
0245   {
0246     if (Verbosity() > 1) { std::cout<<"nmicromegasstates "<<count_clusters<TrkrDefs::micromegasId>(state_keys)<<" , nmicromegasstates cut "<<m_ntpot_states_low_cut<<std::endl; }
0247     return false;
0248   }
0249
0250   return true;
0251 }
0252
0253 int PHTrackPruner::End(PHCompositeNode * /*unused*/)
0254 {
0255   return Fun4AllReturnCodes::EVENT_OK;
0256 }
0257
0258 int PHTrackPruner::GetNodes(PHCompositeNode *topNode)
0259 {
0260   //---------------------------------
0261   // Get additional objects off the Node Tree
0262   //---------------------------------
0263
0264   _svtx_track_map = findNode::getClass<SvtxTrackMap>(topNode, _svtx_track_map_name);
0265   if (!_svtx_track_map)
0266   {
0267     cerr << PHWHERE << " ERROR: Can't find " << _svtx_track_map_name.c_str()  << endl;
0268     return Fun4AllReturnCodes::ABORTEVENT;
0269   }
0270
0271   _si_seed_map = findNode::getClass<TrackSeedContainer>(topNode, _si_seed_map_name);
0272   if (!_si_seed_map)
0273   {
0274     cerr << PHWHERE << " ERROR: Can't find " << _si_seed_map_name.c_str()  << endl;
0275     return Fun4AllReturnCodes::ABORTEVENT;
0276   }
0277
0278   _tpc_seed_map = findNode::getClass<TrackSeedContainer>(topNode, _tpc_seed_map_name);
0279   if (!_tpc_seed_map)
0280   {
0281     cerr << PHWHERE << " ERROR: Can't find " << _tpc_seed_map_name.c_str() << endl;
0282     return Fun4AllReturnCodes::ABORTEVENT;
0283   }
0284
0285   _svtx_seed_map = findNode::getClass<TrackSeedContainer>(topNode, _svtx_seed_map_name);
0286   if (!_svtx_seed_map)
0287   {
0288     cerr << PHWHERE << " ERROR: Can't find " << _svtx_seed_map_name.c_str() << endl;
0289     return Fun4AllReturnCodes::ABORTEVENT;
0290   }
0291
0292   _pruned_svtx_seed_map = findNode::getClass<TrackSeedContainer>(topNode, _pruned_svtx_seed_map_name);
0293   if (!_pruned_svtx_seed_map)
0294   {
0295     std::cout << "Creating node " << _pruned_svtx_seed_map_name.c_str() << std::endl;
0296     /// Get the DST Node
0297     PHNodeIterator iter(topNode);
0298     PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0299
0300     /// Check that it is there
0301     if (!dstNode)
0302     {
0303       std::cerr << "DST Node missing, quitting" << std::endl;
0304       throw std::runtime_error("failed to find DST node in PHTrackPruner::GetNodes");
0305     }
0306
0307     /// Get the tracking subnode
0308     PHNodeIterator dstIter(dstNode);
0309     PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0310
0311     /// Check that it is there
0312     if (!svtxNode)
0313     {
0314       svtxNode = new PHCompositeNode("SVTX");
0315       dstNode->addNode(svtxNode);
0316     }
0317
0318     _pruned_svtx_seed_map = new TrackSeedContainer_v1();
0319     PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(_pruned_svtx_seed_map, _pruned_svtx_seed_map_name, "PHObject");
0320     svtxNode->addNode(node);
0321   }
0322
0323   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0324   if (!_cluster_map)
0325   {
0326     std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
0327     return Fun4AllReturnCodes::ABORTEVENT;
0328   }
0329
0330   _tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0331   if (!_tGeometry)
0332   {
0333     std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
0334     return Fun4AllReturnCodes::ABORTEVENT;
0335   }
0336
0337   return Fun4AllReturnCodes::EVENT_OK;
0338 }
0339
0340 short int PHTrackPruner::findCrossingGeometrically(unsigned int tpcid, unsigned int si_id)
0341 {
0342   // loop over all matches and check for ones with no INTT clusters in the silicon seed
0343   TrackSeed *si_track = _si_seed_map->get(si_id);
0344   const short int crossing = si_track->get_crossing();
0345   const double si_z = TrackSeedHelper::get_z(si_track);
0346
0347   TrackSeed *tpc_track = _tpc_seed_map->get(tpcid);
0348   const double tpc_z = TrackSeedHelper::get_z(tpc_track);
0349
0350   // this is an initial estimate of the bunch crossing based on the z-mismatch of the seeds for this track
0351   short int crossing_estimate = (short int) getBunchCrossing(tpcid, tpc_z - si_z);
0352
0353   if (Verbosity() > 1)
0354   {
0355     std::cout << "findCrossing: "
0356               << " tpcid " << tpcid << " si_id " << si_id << " tpc_z " << tpc_z << " si_z " << si_z << " dz " << tpc_z - si_z
0357               << " INTT crossing " << crossing << " crossing_estimate " << crossing_estimate << std::endl;
0358   }
0359
0360   return crossing_estimate;
0361 }
0362
0363 double PHTrackPruner::getBunchCrossing(unsigned int trid, double z_mismatch)
0364 {
0365   const double vdrift = _tGeometry->get_drift_velocity();  // cm/ns
0366   const double z_bunch_separation = sphenix_constants::time_between_crossings * vdrift; // cm
0367
0368   // The sign of z_mismatch will depend on which side of the TPC the tracklet is in
0369   TrackSeed *track = _tpc_seed_map->get(trid);
0370
0371   // crossing
0372   double crossings = z_mismatch / z_bunch_separation;
0373
0374   // Check the TPC side for the first cluster in the track
0375   unsigned int side = 10;
0376   std::set<short int> side_set;
0377   for (TrackSeed::ConstClusterKeyIter iter = track->begin_cluster_keys();
0378        iter != track->end_cluster_keys();
0379        ++iter)
0380   {
0381     TrkrDefs::cluskey cluster_key = *iter;
0382     unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
0383     if (trkrid == TrkrDefs::tpcId)
0384     {
0385       side = TpcDefs::getSide(cluster_key);
0386       side_set.insert(side);
0387     }
0388   }
0389
0390   if (side == 10)
0391   {
0392     return SHRT_MAX;
0393   }
0394
0395   if (side_set.size() == 2 && Verbosity() > 1)
0396   {
0397     std::cout << "     WARNING: tpc seed " << trid << " changed TPC sides, "
0398               << "  final side " << side << std::endl;
0399   }
0400
0401   // if side = 1 (north, +ve z side), a positive t0 will make the cluster late relative to true z, so it will look like z is less positive
0402   // so a negative z mismatch for side 1 means a positive t0, and positive crossing, so reverse the sign for side 1
0403   if (side == 1)
0404   {
0405     crossings *= -1.0;
0406   }
0407
0408   if (Verbosity() > 1)
0409   {
0410     std::cout << "  gettrackid " << trid << " side " << side << " z_mismatch " << z_mismatch << " crossings " << crossings << std::endl;
0411   }
0412
0413   return crossings;
0414 }