packages/trackreco/PHSiliconTpcTrackMatching.cc

0001 #include "PHSiliconTpcTrackMatching.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 using namespace std;
0043
0044 //____________________________________________________________________________..
0045 PHSiliconTpcTrackMatching::PHSiliconTpcTrackMatching(const std::string &name)
0046   : SubsysReco(name)
0047   , PHParameterInterface(name)
0048 {
0049   InitializeParameters();
0050 }
0051
0052 //____________________________________________________________________________..
0053 PHSiliconTpcTrackMatching::~PHSiliconTpcTrackMatching() = default;
0054
0055 //____________________________________________________________________________..
0056 int PHSiliconTpcTrackMatching::InitRun(PHCompositeNode *topNode)
0057 {
0058   UpdateParametersWithMacro();
0059   if(_test_windows)
0060   {
0061   _file = new TFile(_file_name.c_str(), "RECREATE");
0062   _tree = new TNtuple("track_match", "track_match",
0063                       "event:sicrossing:siq:siphi:sieta:six:siy:siz:sipx:sipy:sipz:tpcq:tpcphi:tpceta:tpcx:tpcy:tpcz:tpcpx:tpcpy:tpcpz:tpcid:siid");
0064   }
0065   // put these in the output file
0066   cout << PHWHERE << " Search windows: phi " << _phi_search_win << " eta "
0067        << _eta_search_win << " _pp_mode " << _pp_mode << " _use_intt_crossing " << _use_intt_crossing << endl;
0068
0069   int ret = GetNodes(topNode);
0070   if (ret != Fun4AllReturnCodes::EVENT_OK)
0071   {
0072     return ret;
0073   }
0074   std::istringstream stringline(m_fieldMap);
0075   stringline >> fieldstrength;
0076
0077   // initialize the WindowMatchers
0078   window_dx.init_bools("dx", _print_windows || Verbosity()   >0);
0079   window_dy.init_bools("dy", _print_windows || Verbosity()   >0);
0080   window_dz.init_bools("dz", _print_windows || Verbosity()   >0);
0081   window_dphi.init_bools("dphi", _print_windows || Verbosity() >0);
0082   window_deta.init_bools("deta", _print_windows || Verbosity() >0);
0083
0084
0085   return ret;
0086 }
0087
0088 //_____________________________________________________________________
0089 void PHSiliconTpcTrackMatching::SetDefaultParameters()
0090 {
0091   // Data on gasses @20 C and 760 Torr from the following source:
0092   // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
0093   // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations:
0094   // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html
0095
0096   return;
0097 }
0098
0099 std::string PHSiliconTpcTrackMatching::WindowMatcher::print_fn(const Arr3D& dat) {
0100   std::ostringstream os;
0101   if (dat[1]==0.) {
0102     os << dat[0];
0103   } else {
0104     os << dat[0] << (dat[1]>0 ? "+" : "") << dat[1] <<"*exp("<<dat[2]<<"/pT)";
0105   }
0106   return os.str();
0107 }
0108
0109 void PHSiliconTpcTrackMatching::WindowMatcher::init_bools(const std::string& tag, const bool print) {
0110   // set values for positive tracks
0111   fabs_max_posQ = (posLo[0]==100.);
0112   posLo_b0 = (posLo[1]==0.);
0113   posHi_b0 = (posHi[1]==0.);
0114   // if no values for negative tracks, copy over from positive tracks
0115   if (negHi[0]==100.) {
0116     negLo = posLo;
0117     negHi = posHi;
0118     fabs_max_negQ = fabs_max_posQ;
0119     negLo_b0 = posLo_b0;
0120     negHi_b0 = posHi_b0;
0121     min_pt_negQ = min_pt_posQ;
0122   } else {
0123     fabs_max_negQ = (negLo[0]==100.);
0124     negLo_b0 = (negLo[1]==0.);
0125     negHi_b0 = (negHi[1]==0.);
0126   }
0127   if (print) {
0128     std::cout << " Track matching window, " << tag << ":" << std::endl;
0129
0130     if (posHi==negHi && posLo == negLo) {
0131       std::cout << "  all tracks: ";
0132     } else {
0133       std::cout << "   +Q tracks: ";
0134     }
0135     if (posLo[0]==100) {
0136       std::cout << "  |" << tag <<"| < " << print_fn(posHi) << std::endl;
0137     } else {
0138       std::cout << print_fn(posLo) <<" < " << tag << " < " << print_fn(posHi) << std::endl;
0139     }
0140
0141     if (posHi != negHi || posLo != negLo) {
0142       std::cout << "   -Q tracks: ";
0143       if (negLo[0]==100) {
0144         std::cout << "  |" << tag <<"| < " << print_fn(negHi) << std::endl;
0145       } else {
0146         std::cout << print_fn(negLo) <<" < " << tag << " < " << print_fn(negHi) << std::endl;
0147       }
0148     }
0149   }
0150
0151 }
0152
0153 bool PHSiliconTpcTrackMatching::WindowMatcher::in_window
0154 (const bool posQ, const double tpc_pt, const double tpc_X, const double si_X)
0155 {
0156   const auto delta = tpc_X-si_X;
0157
0158   if (posQ) {
0159     double pt = (tpc_pt<min_pt_posQ) ? min_pt_posQ : tpc_pt;
0160     if (fabs_max_posQ) {
0161       return fabs(delta) < fn_exp(posHi, posHi_b0, pt);
0162     } else {
0163       return (delta > fn_exp(posLo, posLo_b0, pt)
0164            && delta < fn_exp(posHi, posHi_b0, pt));
0165     }
0166   } else {
0167     double pt = (tpc_pt<min_pt_negQ) ? min_pt_negQ : tpc_pt;
0168     if (fabs_max_negQ) {
0169       return fabs(delta) < fn_exp(negHi, negHi_b0, pt);
0170     } else {
0171       return (delta > fn_exp(negLo, negLo_b0, pt)
0172            && delta < fn_exp(negHi, negHi_b0, pt));
0173     }
0174   }
0175 }
0176
0177 //____________________________________________________________________________..
0178 int PHSiliconTpcTrackMatching::process_event(PHCompositeNode * /*unused*/)
0179 {
0180   if(Verbosity() > 2)
0181   {
0182     std::cout << " Warning: PHSiliconTpcTrackMatching "
0183       << ( _zero_field ? "zero field is ON" : " zero field is OFF") << std::endl;
0184   }
0185   // _track_map contains the TPC seed track stubs
0186   // _track_map_silicon contains the silicon seed track stubs
0187   // _svtx_seed_map contains the combined silicon and tpc track seeds
0188
0189     // in case these objects are in the input file, we clear the nodes and replace them
0190   _svtx_seed_map->Reset();
0191
0192   if (Verbosity() > 0)
0193   {
0194     cout << PHWHERE << " TPC track map size " << _track_map->size() << " Silicon track map size " << _track_map_silicon->size() << endl;
0195   }
0196
0197   if (_track_map->size() == 0)
0198   {
0199     return Fun4AllReturnCodes::EVENT_OK;
0200   }
0201
0202   // loop over the silicon seeds and add the crossing to them
0203   for (unsigned int trackid = 0; trackid != _track_map_silicon->size(); ++trackid)
0204   {
0205     _tracklet_si = _track_map_silicon->get(trackid);
0206     if (!_tracklet_si)
0207     {
0208       continue;
0209     }
0210     auto crossing = _tracklet_si->get_crossing();
0211     if (Verbosity() > 8)
0212     {
0213       std::cout << " silicon stub: " << trackid << " eta " << _tracklet_si->get_eta()
0214         << " pt " << _tracklet_si->get_pt() << " si z " << TrackSeedHelper::get_z(_tracklet_si)
0215         << " crossing " << crossing << std::endl;
0216     }
0217
0218     if (Verbosity() > 1)
0219     {
0220       cout << " Si track " << trackid << " crossing " << crossing << endl;
0221     }
0222   }
0223
0224   // Find all matches of tpc and si tracklets in eta and phi, x and y
0225   std::multimap<unsigned int, unsigned int> tpc_matches;
0226   std::set<unsigned int> tpc_matched_set;
0227   std::set<unsigned int> tpc_unmatched_set;
0228   findEtaPhiMatches(tpc_matched_set, tpc_unmatched_set, tpc_matches);
0229
0230   // check z matching for all matches of tpc and si
0231   // for _pp_mode=false, assume zero crossings for all tracks
0232   // for _pp_mode=true, correct tpc seed z according to crossing number, do nothing if crossing number is not set
0233   // remove matches from tpc_matches if z matching is not satisfied
0234   std::multimap<unsigned int, unsigned int> bad_map;
0235   checkZMatches(tpc_matches, bad_map);
0236
0237   // update tpc_matched_set and tpc_unmatched_set
0238   tpc_matched_set.clear();
0239   for (const auto& [key, _] : tpc_matches)
0240   {
0241     tpc_matched_set.insert(key);
0242   }
0243
0244   for (const auto& [key, _] : bad_map)
0245   {
0246     if (!tpc_matched_set.count(key))
0247     {
0248         tpc_unmatched_set.insert(key);
0249     }
0250   }
0251
0252   // We have a complete list of all eta/phi matched tracks in the map "tpc_matches"
0253   // make the combined track seeds from tpc_matches
0254   for (auto [tpcid, si_id] : tpc_matches)
0255   {
0256     auto svtxseed = std::make_unique<SvtxTrackSeed_v2>();
0257     svtxseed->set_silicon_seed_index(si_id);
0258     svtxseed->set_tpc_seed_index(tpcid);
0259     // In pp mode, if a matched track does not have INTT clusters we have to find the crossing geometrically
0260     // Record the geometrically estimated crossing in the track seeds for later use if needed
0261     short int crossing_estimate = findCrossingGeometrically(tpcid, si_id);
0262     svtxseed->set_crossing_estimate(crossing_estimate);
0263     _svtx_seed_map->insert(svtxseed.get());
0264
0265     if (Verbosity() > 1)
0266     {
0267       std::cout << "  combined seed id " << _svtx_seed_map->size() - 1 << " si id " << si_id << " tpc id " << tpcid << " crossing estimate " << crossing_estimate << std::endl;
0268     }
0269   }
0270
0271   // Also make the unmatched TPC seeds into SvtxTrackSeeds
0272   for (auto tpcid : tpc_unmatched_set)
0273   {
0274     auto svtxseed = std::make_unique<SvtxTrackSeed_v2>();
0275     svtxseed->set_tpc_seed_index(tpcid);
0276     _svtx_seed_map->insert(svtxseed.get());
0277
0278     if (Verbosity() > 1)
0279     {
0280       std::cout << "  converted unmatched TPC seed id " << _svtx_seed_map->size() - 1 << " tpc id " << tpcid << std::endl;
0281     }
0282   }
0283
0284   if (Verbosity() > 0)
0285   {
0286     std::cout << "final svtx seed map size " << _svtx_seed_map->size() << std::endl;
0287   }
0288
0289   if (Verbosity() > 1)
0290   {
0291     for (const auto &seed : *_svtx_seed_map)
0292     {
0293       seed->identify();
0294       std::cout << std::endl;
0295     }
0296
0297     cout << "PHSiliconTpcTrackMatching::process_event(PHCompositeNode *topNode) Leaving process_event" << endl;
0298   }
0299   m_event++;
0300   return Fun4AllReturnCodes::EVENT_OK;
0301 }
0302
0303 short int PHSiliconTpcTrackMatching::findCrossingGeometrically(unsigned int tpcid, unsigned int si_id)
0304 {
0305   // loop over all matches and check for ones with no INTT clusters in the silicon seed
0306   TrackSeed *si_track = _track_map_silicon->get(si_id);
0307   const short int crossing = si_track->get_crossing();
0308   const double si_z = TrackSeedHelper::get_z(si_track);
0309
0310   TrackSeed *tpc_track = _track_map->get(tpcid);
0311   const double tpc_z = TrackSeedHelper::get_z(tpc_track);
0312
0313   // this is an initial estimate of the bunch crossing based on the z-mismatch of the seeds for this track
0314   short int crossing_estimate = (short int) getBunchCrossing(tpcid, tpc_z - si_z);
0315
0316   if (Verbosity() > 1)
0317   {
0318     std::cout << "findCrossing: "
0319               << " tpcid " << tpcid << " si_id " << si_id << " tpc_z " << tpc_z << " si_z " << si_z << " dz " << tpc_z - si_z
0320               << " INTT crossing " << crossing << " crossing_estimate " << crossing_estimate << std::endl;
0321   }
0322
0323   return crossing_estimate;
0324 }
0325
0326 double PHSiliconTpcTrackMatching::getBunchCrossing(unsigned int trid, double z_mismatch)
0327 {
0328   const double vdrift = _tGeometry->get_drift_velocity();  // cm/ns
0329   const double z_bunch_separation = sphenix_constants::time_between_crossings * vdrift; // cm
0330
0331   // The sign of z_mismatch will depend on which side of the TPC the tracklet is in
0332   TrackSeed *track = _track_map->get(trid);
0333
0334   // crossing
0335   double crossings = z_mismatch / z_bunch_separation;
0336
0337   // Check the TPC side for the first cluster in the track
0338   unsigned int side = 10;
0339   std::set<short int> side_set;
0340   for (TrackSeed::ConstClusterKeyIter iter = track->begin_cluster_keys();
0341        iter != track->end_cluster_keys();
0342        ++iter)
0343   {
0344     TrkrDefs::cluskey cluster_key = *iter;
0345     unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
0346     if (trkrid == TrkrDefs::tpcId)
0347     {
0348       side = TpcDefs::getSide(cluster_key);
0349       side_set.insert(side);
0350     }
0351   }
0352
0353   if (side == 10)
0354   {
0355     return SHRT_MAX;
0356   }
0357
0358   if (side_set.size() == 2 && Verbosity() > 1)
0359   {
0360     std::cout << "     WARNING: tpc seed " << trid << " changed TPC sides, "
0361               << "  final side " << side << std::endl;
0362   }
0363
0364   // 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
0365   // so a negative z mismatch for side 1 means a positive t0, and positive crossing, so reverse the sign for side 1
0366   if (side == 1)
0367   {
0368     crossings *= -1.0;
0369   }
0370
0371   if (Verbosity() > 1)
0372   {
0373     std::cout << "  gettrackid " << trid << " side " << side << " z_mismatch " << z_mismatch << " crossings " << crossings << std::endl;
0374   }
0375
0376   return crossings;
0377 }
0378
0379 int PHSiliconTpcTrackMatching::End(PHCompositeNode * /*unused*/)
0380 {
0381   if(_test_windows)
0382   {
0383   _file->cd();
0384   _tree->Write();
0385   _file->Close();
0386   }
0387   return Fun4AllReturnCodes::EVENT_OK;
0388 }
0389
0390 int PHSiliconTpcTrackMatching::GetNodes(PHCompositeNode *topNode)
0391 {
0392   //---------------------------------
0393   // Get additional objects off the Node Tree
0394   //---------------------------------
0395
0396   _cluster_crossing_map = findNode::getClass<TrkrClusterCrossingAssoc>(topNode, "TRKR_CLUSTERCROSSINGASSOC");
0397   if (!_cluster_crossing_map)
0398   {
0399     //cerr << PHWHERE << " ERROR: Can't find TRKR_CLUSTERCROSSINGASSOC " << endl;
0400     // return Fun4AllReturnCodes::ABORTEVENT;
0401   }
0402
0403   _track_map_silicon = findNode::getClass<TrackSeedContainer>(topNode, _silicon_track_map_name);
0404   if (!_track_map_silicon)
0405   {
0406     cerr << PHWHERE << " ERROR: Can't find SiliconTrackSeedContainer " << endl;
0407     return Fun4AllReturnCodes::ABORTEVENT;
0408   }
0409
0410   _track_map = findNode::getClass<TrackSeedContainer>(topNode, _track_map_name);
0411   if (!_track_map)
0412   {
0413     cerr << PHWHERE << " ERROR: Can't find " << _track_map_name.c_str() << endl;
0414     return Fun4AllReturnCodes::ABORTEVENT;
0415   }
0416
0417   _svtx_seed_map = findNode::getClass<TrackSeedContainer>(topNode, "SvtxTrackSeedContainer");
0418   if (!_svtx_seed_map)
0419   {
0420     std::cout << "Creating node SvtxTrackSeedContainer" << std::endl;
0421     /// Get the DST Node
0422     PHNodeIterator iter(topNode);
0423     PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0424
0425     /// Check that it is there
0426     if (!dstNode)
0427     {
0428       std::cerr << "DST Node missing, quitting" << std::endl;
0429       throw std::runtime_error("failed to find DST node in PHActsSourceLinks::createNodes");
0430     }
0431
0432     /// Get the tracking subnode
0433     PHNodeIterator dstIter(dstNode);
0434     PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0435
0436     /// Check that it is there
0437     if (!svtxNode)
0438     {
0439       svtxNode = new PHCompositeNode("SVTX");
0440       dstNode->addNode(svtxNode);
0441     }
0442
0443     _svtx_seed_map = new TrackSeedContainer_v1();
0444     PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(_svtx_seed_map, "SvtxTrackSeedContainer", "PHObject");
0445     svtxNode->addNode(node);
0446   }
0447
0448   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0449   if (!_cluster_map)
0450   {
0451     std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
0452     return Fun4AllReturnCodes::ABORTEVENT;
0453   }
0454
0455   _tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0456   if (!_tGeometry)
0457   {
0458     std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
0459     return Fun4AllReturnCodes::ABORTEVENT;
0460   }
0461
0462   return Fun4AllReturnCodes::EVENT_OK;
0463 }
0464
0465 void PHSiliconTpcTrackMatching::findEtaPhiMatches(
0466     std::set<unsigned int> &tpc_matched_set,
0467     std::set<unsigned int> &tpc_unmatched_set,
0468     std::multimap<unsigned int, unsigned int> &tpc_matches)
0469 {
0470   // loop over the TPC track seeds
0471   for (unsigned int phtrk_iter = 0;
0472        phtrk_iter < _track_map->size();
0473        ++phtrk_iter)
0474   {
0475     _tracklet_tpc = _track_map->get(phtrk_iter);
0476     if (!_tracklet_tpc)
0477     {
0478       continue;
0479     }
0480
0481     unsigned int tpcid = phtrk_iter;
0482     if (Verbosity() > 1)
0483     {
0484       std::cout
0485           << __LINE__
0486           << ": Processing seed itrack: " << tpcid
0487           << ": nhits: " << _tracklet_tpc->size_cluster_keys()
0488           << ": Total tracks: " << _track_map->size()
0489           << ": phi: " << _tracklet_tpc->get_phi()
0490           << endl;
0491     }
0492
0493     double tpc_phi, tpc_eta, tpc_pt;
0494     float tpc_px, tpc_py, tpc_pz;
0495     int tpc_q;
0496     Acts::Vector3 tpc_pos;
0497     if (_zero_field) {
0498       auto cluster_list = getTrackletClusterList(_tracklet_tpc);
0499
0500       Acts::Vector3  mom;
0501       bool ok_track;
0502
0503       std::tie(ok_track, tpc_phi, tpc_eta, tpc_pt, tpc_pos, mom) =
0504         TrackFitUtils::zero_field_track_params(_tGeometry, _cluster_map, cluster_list);
0505       if (!ok_track) { continue; }
0506       tpc_px = mom.x();
0507       tpc_py = mom.y();
0508       tpc_pz = mom.z();
0509       tpc_q = -100;
0510     } else {
0511       tpc_phi = _tracklet_tpc->get_phi();
0512       tpc_eta = _tracklet_tpc->get_eta();
0513       tpc_pt = fabs(1. / _tracklet_tpc->get_qOverR()) * (0.3 / 100.) * fieldstrength;
0514
0515       tpc_pos = TrackSeedHelper::get_xyz(_tracklet_tpc);
0516
0517       tpc_px = _tracklet_tpc->get_px();
0518       tpc_py = _tracklet_tpc->get_py();
0519       tpc_pz = _tracklet_tpc->get_pz();
0520
0521       tpc_q = _tracklet_tpc->get_charge();
0522     }
0523
0524     bool is_posQ = (tpc_q>0.);
0525
0526     if (Verbosity() > 8)
0527     {
0528       std::cout << " tpc stub: " << tpcid << " eta " << tpc_eta << " phi " << tpc_phi << " pt " << tpc_pt << " tpc z " << TrackSeedHelper::get_z(_tracklet_tpc) << std::endl;
0529     }
0530
0531     if (Verbosity() > 3)
0532     {
0533       cout << "TPC tracklet:" << endl;
0534       _tracklet_tpc->identify();
0535     }
0536
0537
0538     bool matched = false;
0539
0540     // Now search the silicon track list for a match in eta and phi
0541     for (unsigned int phtrk_iter_si = 0;
0542          phtrk_iter_si < _track_map_silicon->size();
0543          ++phtrk_iter_si)
0544     {
0545       _tracklet_si = _track_map_silicon->get(phtrk_iter_si);
0546       if (!_tracklet_si)
0547       {
0548
0549         continue;
0550       }
0551
0552       double si_phi, si_eta, si_pt;
0553       float si_px, si_py, si_pz;
0554       int si_q;
0555       Acts::Vector3 si_pos;
0556       if (_zero_field) {
0557       auto cluster_list = getTrackletClusterList(_tracklet_si);
0558
0559       Acts::Vector3  mom;
0560       bool ok_track;
0561
0562       std::tie(ok_track, si_phi, si_eta, si_pt, si_pos, mom) =
0563         TrackFitUtils::zero_field_track_params(_tGeometry, _cluster_map, cluster_list);
0564       if (!ok_track) { continue; }
0565         si_px = mom.x();
0566         si_py = mom.y();
0567         si_pz = mom.z();
0568         si_q = -100;
0569       } else {
0570         si_eta = _tracklet_si->get_eta();
0571         si_phi = _tracklet_si->get_phi();
0572
0573         si_pos = TrackSeedHelper::get_xyz(_tracklet_si);
0574         si_px = _tracklet_si->get_px();
0575         si_py = _tracklet_si->get_py();
0576         si_pz = _tracklet_si->get_pz();
0577         si_q = _tracklet_si->get_charge();
0578       }
0579       int si_crossing = _tracklet_si->get_crossing();
0580       unsigned int siid = phtrk_iter_si;
0581
0582       if(_test_windows)
0583       {
0584         float data[] = {
0585           (float) m_event, (float) si_crossing,
0586           (float) si_q, (float) si_phi, (float) si_eta, (float) si_pos.x(), (float) si_pos.y(), (float) si_pos.z(), (float) si_px, (float) si_py, (float) si_pz,
0587           (float) tpc_q, (float) tpc_phi, (float) tpc_eta, (float) tpc_pos.x(), (float) tpc_pos.y(), (float) tpc_pos.z(), (float) tpc_px, (float) tpc_py, (float) tpc_pz,
0588           (float) tpcid, (float) siid
0589     };
0590         _tree->Fill(data);
0591       }
0592
0593       bool eta_match = false;
0594       if (window_deta.in_window(is_posQ, tpc_pt, tpc_eta, si_eta))
0595       {
0596         eta_match = true;
0597       }
0598       else if (fabs(tpc_eta-si_eta) < _deltaeta_min)
0599       {
0600     eta_match = true;
0601       }
0602       if (!eta_match)
0603       {
0604         continue;
0605       }
0606
0607       bool position_match = false;
0608       if (window_dx.in_window(is_posQ, tpc_pt, tpc_pos.x(), si_pos.x())
0609        && window_dy.in_window(is_posQ, tpc_pt, tpc_pos.y(), si_pos.y()))
0610       {
0611         position_match = true;
0612       }
0613       if (!position_match)
0614       {
0615         continue;
0616       }
0617
0618       bool phi_match = false;
0619       if (window_dphi.in_window(is_posQ, tpc_pt, tpc_phi, si_phi))
0620       {
0621         phi_match = true;
0622         // if phi fails, account for case where |tpc_phi-si_phi|>PI
0623       } else if (fabs(tpc_phi-si_phi)>M_PI) {
0624         auto tpc_phi_wrap = tpc_phi;
0625         if ((tpc_phi_wrap - si_phi) > M_PI) {
0626           tpc_phi_wrap -= 2*M_PI;
0627         } else {
0628           tpc_phi_wrap += 2*M_PI;
0629         }
0630         phi_match = window_dphi.in_window(is_posQ, tpc_pt, tpc_phi_wrap, si_phi);
0631       }
0632       if (!phi_match)
0633       {
0634         continue;
0635       }
0636
0637       if (Verbosity() > 3)
0638       {
0639         cout << " testing for a match for TPC track " << tpcid << " with pT " << _tracklet_tpc->get_pt()
0640              << " and eta " << _tracklet_tpc->get_eta() << " with Si track " << siid << " with crossing " << _tracklet_si->get_crossing() << endl;
0641         cout << " tpc_phi " << tpc_phi << " si_phi " << si_phi << " dphi " << tpc_phi - si_phi << " phi search " << _phi_search_win << " tpc_eta " << tpc_eta
0642              << " si_eta " << si_eta << " deta " << tpc_eta - si_eta << " eta search " << _eta_search_win  << endl;
0643         std::cout << "      tpc x " << tpc_pos.x() << " si x " << si_pos.x() << " tpc y " << tpc_pos.y() << " si y " << si_pos.y() << " tpc_z " << tpc_pos.z() << " si z " << si_pos.z() << std::endl;
0644         std::cout << "      x search " << _x_search_win  << " y search " << _y_search_win << " z search " << _z_search_win << std::endl;
0645       }
0646
0647       // got a match, add to the list
0648       // These stubs are matched in eta, phi, x and y already
0649       matched = true;
0650       tpc_matches.insert(std::make_pair(tpcid, siid));
0651       tpc_matched_set.insert(tpcid);
0652
0653       if (Verbosity() > 1)
0654       {
0655         cout << " found a match for TPC track " << tpcid << " with Si track " << siid << endl;
0656         cout << "          tpc_phi " << tpc_phi << " si_phi " << si_phi << " phi_match " << phi_match
0657              << " tpc_eta " << tpc_eta << " si_eta " << si_eta << " eta_match " << eta_match << endl;
0658         std::cout << "      tpc x " << tpc_pos.x() << " si x " << si_pos.x() << " tpc y " << tpc_pos.y() << " si y " << si_pos.y() << " tpc_z " << tpc_pos.z() << " si z " << si_pos.z() << std::endl;
0659       }
0660
0661       // temporary!
0662       if (_test_windows && Verbosity() > 1)
0663       {
0664         cout << " Try_silicon: crossing" << si_crossing <<  "  pt " << tpc_pt << " tpc_phi " << tpc_phi << " si_phi " << si_phi << " dphi " << tpc_phi - si_phi <<  "   si_q" << si_q << "   tpc_q" << tpc_q
0665              << " tpc_eta " << tpc_eta << " si_eta " << si_eta << " deta " << tpc_eta - si_eta << " tpc_x " << tpc_pos.x() << " tpc_y " << tpc_pos.y() << " tpc_z " << tpc_pos.z()
0666              << " dx " << tpc_pos.x() - si_pos.x() << " dy " << tpc_pos.y() - si_pos.y() << " dz " << tpc_pos.z() - si_pos.z()
0667              << endl;
0668       }
0669
0670     }
0671     // if no match found, keep tpc seed for fitting
0672     if (!matched)
0673     {
0674       if (Verbosity() > 1)
0675       {
0676         cout << "inserted unmatched tpc seed " << tpcid << endl;
0677       }
0678       tpc_unmatched_set.insert(tpcid);
0679     }
0680   }
0681
0682   return;
0683 }
0684 void PHSiliconTpcTrackMatching::checkZMatches(
0685     std::multimap<unsigned int, unsigned int> &tpc_matches,
0686     std::multimap<unsigned int, unsigned int> &bad_map)
0687 {
0688   // for _pp_mode=false, assume zero crossings for all track matches
0689   // for _pp_mode=true, do crossing correction on track position z according to side and vdrift
0690   // z matching criteria follows window_z
0691   // there is a dz threshold cut to avoid window_z blow up at low pT
0692
0693   float vdrift = _tGeometry->get_drift_velocity();
0694
0695   for (auto [tpcid, si_id] : tpc_matches)
0696   {
0697     TrackSeed *tpc_track = _track_map->get(tpcid);
0698     TrackSeed *si_track = _track_map_silicon->get(si_id);
0699
0700     short int crossing = si_track->get_crossing();
0701     float tpc_pt, tpc_z, si_z;
0702     int tpc_q;
0703     if (_zero_field) {
0704       auto cluster_list_tpc = getTrackletClusterList(tpc_track);
0705       auto cluster_list_si = getTrackletClusterList(si_track);
0706
0707       tpc_pt = std::get<3>(TrackFitUtils::zero_field_track_params(_tGeometry, _cluster_map, cluster_list_tpc));
0708       tpc_z = std::get<4>(TrackFitUtils::zero_field_track_params(_tGeometry, _cluster_map, cluster_list_tpc)).z();
0709       tpc_q = -100;
0710
0711       si_z = std::get<4>(TrackFitUtils::zero_field_track_params(_tGeometry, _cluster_map, cluster_list_si)).z();
0712     } else {
0713       tpc_pt = fabs(1. / _tracklet_tpc->get_qOverR()) * (0.3 / 100.) * fieldstrength;
0714       tpc_z = TrackSeedHelper::get_z(tpc_track);
0715       tpc_q = _tracklet_tpc->get_charge();
0716       si_z = TrackSeedHelper::get_z(si_track);
0717     }
0718
0719     // get TPC side from one of the TPC clusters
0720     std::vector<TrkrDefs::cluskey> temp_clusters = getTrackletClusterList(tpc_track);
0721     if(temp_clusters.size() == 0) { continue; }
0722     unsigned int this_side = TpcDefs::getSide(temp_clusters[0]);
0723
0724     bool is_posQ = (tpc_q>0.);
0725
0726     float z_mismatch = tpc_z - si_z;
0727     float tpc_z_corrected = _clusterCrossingCorrection.correctZ(tpc_z, this_side, crossing);
0728     float z_mismatch_corrected = tpc_z_corrected - si_z;
0729
0730     bool z_match = false;
0731     if (_pp_mode)
0732     {
0733       if (crossing == SHRT_MAX)
0734       {
0735         if (Verbosity() > 2)
0736         {
0737           std::cout << " drop si_track " << si_id << " with eta " << si_track->get_eta() << " and z " << TrackSeedHelper::get_z(si_track) << " because crossing is undefined " << std::endl;
0738         }
0739         continue;
0740       }
0741
0742       if (window_dz.in_window(is_posQ, tpc_pt, tpc_z_corrected, si_z) && (fabs(z_mismatch_corrected) < _crossing_deltaz_max))
0743       {
0744         z_match = true;
0745       }
0746       else if (fabs(z_mismatch_corrected) < _crossing_deltaz_min)
0747       {
0748     z_match = true;
0749       }
0750     }
0751     else
0752     {
0753       if (window_dz.in_window(is_posQ, tpc_pt, tpc_z, si_z) && (fabs(z_mismatch) < _crossing_deltaz_max))
0754       {
0755         z_match = true;
0756       }
0757       else if (fabs(z_mismatch) < _crossing_deltaz_min)
0758       {
0759     z_match = true;
0760       }
0761     }
0762
0763     if (z_match)
0764     {
0765       if (Verbosity() > 1)
0766       {
0767         std::cout << "  Success:  crossing " << crossing << " tpcid " << tpcid << " si id " << si_id
0768                   << " tpc z " << tpc_z << " si z " << si_z << " z_mismatch " << z_mismatch << "tpc z corrected " << tpc_z_corrected
0769                   << " z_mismatch_corrected " << z_mismatch_corrected << " drift velocity " << vdrift << std::endl;
0770       }
0771     }
0772     else
0773     {
0774       if (Verbosity() > 1)
0775       {
0776         std::cout << "  FAILURE:  crossing " << crossing << " tpcid " << tpcid << " si id " << si_id
0777                   << " tpc z " << tpc_z << " si z " << si_z << " z_mismatch " << z_mismatch << "tpc_z_corrected " << tpc_z_corrected
0778                   << " z_mismatch_corrected " << z_mismatch_corrected << std::endl;
0779       }
0780
0781       bad_map.insert(std::make_pair(tpcid, si_id));
0782     }
0783   }
0784
0785   // remove bad entries from tpc_matches
0786   for (auto [tpcid, si_id] : bad_map)
0787   {
0788     // Have to iterate over tpc_matches and examine each pair to find the one matching bad_map
0789     // this logic works because we call the equal range on vertex_map for every id_pair
0790     // so we only delete one entry per equal range call
0791     auto ret = tpc_matches.equal_range(tpcid);
0792     for (auto it = ret.first; it != ret.second; ++it)
0793     {
0794       if (it->first == tpcid && it->second == si_id)
0795       {
0796         if (Verbosity() > 1)
0797         {
0798           std::cout << "                        erasing tpc_matches entry for tpcid " << tpcid << " si_id " << si_id << std::endl;
0799         }
0800         tpc_matches.erase(it);
0801         break;  // the iterator is no longer valid
0802       }
0803     }
0804   }
0805
0806   return;
0807 }
0808
0809 std::vector<TrkrDefs::cluskey> PHSiliconTpcTrackMatching::getTrackletClusterList(TrackSeed* tracklet)
0810 {
0811   std::vector<TrkrDefs::cluskey> cluskey_vec;
0812   for (auto clusIter = tracklet->begin_cluster_keys();
0813        clusIter != tracklet->end_cluster_keys();
0814        ++clusIter)
0815   {
0816     auto key = *clusIter;
0817     auto cluster = _cluster_map->findCluster(key);
0818     if (!cluster)
0819     {
0820       if(Verbosity() > 1)
0821       {
0822         std::cout << PHWHERE << "Failed to get cluster with key " << key << std::endl;
0823       }
0824       continue;
0825     }
0826
0827     /// Make a safety check for clusters that couldn't be attached to a surface
0828     auto surf = _tGeometry->maps().getSurface(key, cluster);
0829     if (!surf)
0830     {
0831       continue;
0832     }
0833
0834     // drop some bad layers in the TPC completely
0835     unsigned int layer = TrkrDefs::getLayer(key);
0836     if (layer == 7 || layer == 22 || layer == 23 || layer == 38 || layer == 39)
0837     {
0838       continue;
0839     }
0840
0841     cluskey_vec.push_back(key);
0842   }  // end loop over clusters for this track
0843   return cluskey_vec;
0844 }