packages/trackreco/PHTrackSelector.cc

0001 #include "PHTrackSelector.h"
0002
0003 /// Tracking includes
0004
0005 #include <trackbase/TrkrCluster.h>            // for TrkrCluster
0006 #include <trackbase/TrkrDefs.h>               // for cluskey, getLayer, TrkrId
0007 #include <trackbase/TrkrClusterContainer.h>
0008 #include <trackbase/TrkrClusterIterationMapv1.h>
0009 #include <trackbase_historic/SvtxTrack.h>     // for SvtxTrack, SvtxTrack::C...
0010 #include <trackbase_historic/SvtxTrackMap.h>
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/PHIODataNode.h>                         // for PHIODataNode
0016 #include <phool/PHNode.h>                               // for PHNode
0017 #include <phool/PHNodeIterator.h>
0018 #include <phool/PHObject.h>                             // for PHObject
0019 #include <phool/getClass.h>
0020 #include <phool/phool.h>
0021
0022 #include <cmath>                              // for sqrt, fabs, atan2, cos
0023 #include <iostream>                           // for operator<<, basic_ostream
0024 #include <map>                                // for map
0025 #include <set>                                // for _Rb_tree_const_iterator
0026 #include <utility>                            // for pair, make_pair
0027
0028 //____________________________________________________________________________..
0029 PHTrackSelector::PHTrackSelector(const std::string &name)
0030   : SubsysReco(name)
0031   , _track_map_name("SvtxTrackMap")
0032   , _n_iter(1)
0033   , min_tpc_clusters(40)
0034   , min_mvtx_hits(2)
0035   , min_intt_hits(1)
0036   , max_chi2_ndf(30)
0037 {
0038
0039 }
0040
0041 //____________________________________________________________________________..
0042 PHTrackSelector::~PHTrackSelector()
0043 {
0044
0045 }
0046
0047 //____________________________________________________________________________..
0048 int PHTrackSelector::InitRun(PHCompositeNode *topNode)
0049 {
0050   int ret = GetNodes(topNode);
0051   if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
0052
0053   return ret;
0054 }
0055
0056 //____________________________________________________________________________..
0057 int PHTrackSelector::process_event(PHCompositeNode */*topNode*/)
0058 {
0059
0060   if(Verbosity() > 0)
0061     std::cout << PHWHERE << " track map size " << _track_map->size() << std::endl;
0062
0063   if(_n_iter <=1) _iteration_map->Reset();
0064
0065   std::set<unsigned int> track_delete_list;
0066
0067   unsigned int good_track = 0;  // for diagnostic output only
0068   unsigned int bad_track = 0;   // tracks to keep
0069
0070   // make a list of tracks that did not make the keep list
0071   for(auto track_it = _track_map->begin(); track_it != _track_map->end(); ++track_it)
0072     {
0073       auto id = track_it->first;
0074       _track = track_it->second;
0075       bool delete_track = false;
0076       double chi2_ndf = _track->get_chisq() / _track->get_ndf();
0077       unsigned int ntpc = 0;
0078       unsigned int nintt = 0;
0079       unsigned int nmvtx = 0;
0080
0081       for (SvtxTrack::ConstClusterKeyIter iter = _track->begin_cluster_keys();
0082        iter != _track->end_cluster_keys();
0083        ++iter){
0084     TrkrDefs::cluskey cluster_key = *iter;
0085     //  unsigned int layer = TrkrDefs::getLayer(cluster_key);
0086
0087     if(TrkrDefs::getTrkrId(cluster_key) == TrkrDefs::mvtxId ) nmvtx++;
0088     if(TrkrDefs::getTrkrId(cluster_key) == TrkrDefs::inttId ) nintt++;
0089     if(TrkrDefs::getTrkrId(cluster_key) == TrkrDefs::tpcId )  ntpc++;
0090
0091       }
0092
0093       if(chi2_ndf > max_chi2_ndf) delete_track = true;
0094       if(ntpc < min_tpc_clusters) delete_track = true;
0095       if(ntpc < min_tpc_clusters) delete_track = true;
0096       if(nintt < min_intt_hits) delete_track = true;
0097       if(nmvtx < min_mvtx_hits) delete_track = true;
0098
0099       if(delete_track){
0100     track_delete_list.insert(id);
0101     bad_track++;
0102       }
0103       else{
0104     //add clusters to cluster used list
0105     for (SvtxTrack::ConstClusterKeyIter iter = _track->begin_cluster_keys();
0106          iter != _track->end_cluster_keys();
0107          ++iter){
0108       TrkrDefs::cluskey cluster_key = *iter;
0109       _iteration_map->addIteration(cluster_key,_n_iter);
0110     }
0111     good_track++;
0112       }
0113
0114     }
0115   if(Verbosity() > 0)
0116     std::cout << " Number of good tracks is " << good_track << " bad tracks " << bad_track << std::endl;
0117
0118   // delete failed tracks
0119   for(auto it = track_delete_list.begin(); it != track_delete_list.end(); ++it){
0120     if(Verbosity() > 1)
0121       std::cout << " erasing track ID " << *it << std::endl;
0122     _track_map->erase(*it);
0123   }
0124
0125   if(Verbosity() > 0)
0126     std::cout << " track_delete_list size " << track_delete_list.size() << std::endl;
0127
0128   // delete failed tracks
0129   for(auto it = track_delete_list.begin(); it != track_delete_list.end(); ++it)
0130     {
0131       if(Verbosity() > 1)
0132     std::cout << " erasing track ID " << *it << std::endl;
0133       _track_map->erase(*it);
0134     }
0135
0136   if(Verbosity() > 0)
0137     std::cout << "Track map size after choosing best silicon match: " << _track_map->size() << std::endl;
0138
0139   return Fun4AllReturnCodes::EVENT_OK;
0140 }
0141
0142 int PHTrackSelector::End(PHCompositeNode */*topNode*/)
0143 {
0144   return Fun4AllReturnCodes::EVENT_OK;
0145 }
0146
0147 int  PHTrackSelector::GetNodes(PHCompositeNode* topNode)
0148 {
0149   PHNodeIterator iter(topNode);
0150   // Looking for the DST node
0151   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0152   if (!dstNode)
0153   {
0154     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0155     return Fun4AllReturnCodes::ABORTRUN;
0156   }
0157
0158   //  TrackClusterIterationMapv1* _iteration_map;
0159   _track_map = findNode::getClass<SvtxTrackMap>(topNode,  _track_map_name);
0160   if (!_track_map)
0161   {
0162     std::cout << PHWHERE << " ERROR: Can't find SvtxTrackMap: " << std::endl;
0163     return Fun4AllReturnCodes::ABORTEVENT;
0164   }
0165
0166   // Create the Cluster Iteration Map node if required
0167   _iteration_map= findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
0168   if (!_iteration_map)
0169   {
0170     PHNodeIterator dstiter(dstNode);
0171     PHCompositeNode *DetNode =
0172       dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0173     if (!DetNode)
0174     {
0175       DetNode = new PHCompositeNode("TRKR");
0176       dstNode->addNode(DetNode);
0177     }
0178
0179     _iteration_map = new TrkrClusterIterationMapv1;
0180     PHIODataNode<PHObject> *TrkrClusterIterationMapv1 =
0181       new PHIODataNode<PHObject>(_iteration_map, "CLUSTER_ITERATION_MAP", "PHObject");
0182     DetNode->addNode(TrkrClusterIterationMapv1);
0183   }
0184
0185
0186
0187   return Fun4AllReturnCodes::EVENT_OK;
0188 }