27 #if __cplusplus < 201402L
28 #include <boost/make_unique.hpp>
45 , _track_map_name_silicon(
"SvtxSiliconTrackMap")
64 fdphi =
new TF1(
"f1",
"[0] + [1]/x^[2]");
77 fscdphi =
new TF1(
"f2",
"[0] + [1]*x^2");
116 cout <<
" Si track " <<
_tracklet_si->
get_id() <<
" si_phi " << si_phi <<
" si_eta " << si_eta << endl;
124 std::multimap<unsigned int, unsigned int> tpc_matches;
125 std::set<unsigned int> tpc_matched_set;
130 std::multimap<int, std::pair<unsigned int, unsigned int>> crossing_matches;
131 std::map<unsigned int, int> tpc_crossing_map;
132 std::set<int> crossing_set;
140 tagMatchCrossing(tpc_matches, crossing_set, crossing_matches, tpc_crossing_map);
144 std::multimap<double, std::pair<unsigned int, unsigned int>> si_sorted_map;
145 for(
auto ncross : crossing_set)
147 if(
Verbosity() > 1) std::cout <<
" ncross = " << ncross << std::endl;
149 auto ret = crossing_matches. equal_range(ncross);
150 for(
auto it = ret.first;
it != ret.second; ++
it)
153 std::cout <<
" crossing " <<
it->first <<
" tpc_id " <<
it->second.first <<
" si_id " <<
it->second.second
161 si_sorted_map.insert(std::make_pair(z_si, std::make_pair(
it->second.first,
it->second.second)));
166 std::multimap<unsigned int, std::pair<unsigned int, unsigned int>> vertex_map;
167 std::vector<double> vertex_list;
171 std::map<unsigned int, double> vertex_crossings_map;
175 cleanVertexMap( vertex_crossings_map, vertex_map, tpc_crossing_map );
200 cout <<
"PHSiliconTpcTrackMatching::process_event(PHCompositeNode *topNode) Leaving process_event" << endl;
216 _assoc_container = findNode::getClass<AssocInfoContainer>(topNode,
"AssocInfoContainer");
219 cerr <<
PHWHERE <<
" ERROR: Can't find AssocInfoContainer." << endl;
226 cerr <<
PHWHERE <<
" ERROR: Can't find SvtxSiliconTrackMap: " << endl;
233 cerr <<
PHWHERE <<
" ERROR: Can't find SvtxTrackMap " << endl;
240 std::cout <<
"Creating node TpcSeedTrackMap" << std::endl;
249 std::cerr <<
"DST Node missing, quitting" << std::endl;
250 throw std::runtime_error(
"failed to find DST node in PHActsSourceLinks::createNodes");
272 std::cout <<
" Found CORRECTED_TRKR_CLUSTER node " << std::endl;
275 _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
278 std::cout <<
PHWHERE <<
" ERROR: Can't find node TRKR_CLUSTER" << std::endl;
282 _tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
285 std::cout <<
PHWHERE <<
"Error, can't find acts tracking geometry" << std::endl;
289 _surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
292 std::cout <<
PHWHERE <<
"Error, can't find acts surface maps" << std::endl;
302 double vdrift = 8.00;
303 double z_bunch_separation = 0.107 * vdrift;
308 double crossings = z_mismatch / z_bunch_separation;
342 if(clus_z -z_mismatch > 0)
345 if(
Verbosity() > 1) std::cout <<
" trackid " << trid <<
" clus_z " << clus_z <<
" z_mismatch " << z_mismatch <<
" crossings " << crossings << std::endl;
352 if(v.size() == 0)
return NAN;
356 if( (v.size() % 2) == 0)
360 auto m1 = v.begin() + v.size()/2;
361 std::nth_element(v.begin(), m1, v.end());
362 double median1 = v[v.size()/2];
364 auto m2 = v.begin() + v.size()/2 - 1;
365 std::nth_element(v.begin(),
m2, v.end());
366 double median2 = v[v.size()/2 - 1];
368 median = (median1 + median2) / 2.0;
369 if(
Verbosity() > 2) std::cout <<
"The vector size is " << v.size()
370 <<
" element m is " << v.size() / 2 <<
" = " << v[v.size()/2]
371 <<
" element m-1 is " << v.size() / 2 -1 <<
" = " << v[v.size()/2-1]
377 auto m = v.begin() + v.size()/2;
378 std::nth_element(v.begin(),
m, v.end());
379 median = v[v.size()/2];
380 if(
Verbosity() > 2) std::cout <<
"The vector size is " << v.size() <<
" element m is " << v.size() / 2 <<
" = " << v[v.size()/2] << std::endl;
389 for(
auto it = tpc_matches.begin();
it != tpc_matches.end(); ++
it)
391 unsigned int tpcid =
it->first;
393 if(
Verbosity() > 1) std::cout <<
" tpcid " << tpcid <<
" original z " << tpc_track->
get_z() << std::endl;
396 unsigned int si_id =
it->second;
398 if(
Verbosity() > 1) std::cout <<
" si track id " << si_id << std::endl;
406 cout <<
" inserting si cluster key " << si_cluster_key <<
" into existing TPC track " << tpc_track->
get_id() << endl;
420 std::cout <<
" TPC seed track ID " << tpc_track->
get_id() <<
" si track id " << si_track->
get_id()
434 for(
auto it = crossing_matches.begin();
it != crossing_matches.end(); ++
it)
436 unsigned int tpcid =
it->second.first;
438 if(
Verbosity() > 1) std::cout <<
" tpcid " << tpcid <<
" original z " << tpc_track->
get_z() << std::endl;
441 unsigned int si_id =
it->second.second;
443 if(
Verbosity() > 1) std::cout <<
" si track id " << si_id << std::endl;
451 cout <<
" inserting si cluster key " << si_cluster_key <<
" into existing TPC track " << tpc_track->
get_id() << endl;
464 std::cout <<
" TPC seed track ID " << tpc_track->
get_id() <<
" si track id " << si_track->
get_id()
477 for(
auto it = vertex_map.begin();
it != vertex_map.end(); ++
it)
479 unsigned int tpcid =
it->second.first;
481 if(
Verbosity() > 1) std::cout <<
" tpcid " << tpcid <<
" original z " << tpc_track->
get_z() << std::endl;
484 unsigned int si_id =
it->second.second;
486 if(
Verbosity() > 1) std::cout <<
" si track id " << si_id << std::endl;
494 cout <<
" inserting si cluster key " << si_cluster_key <<
" into existing TPC track " << tpc_track->
get_id() << endl;
508 std::cout <<
" TPC seed track ID " << tpc_track->
get_id() <<
" si track id " << si_track->
get_id()
520 std::map<unsigned int, double> &vertex_crossings_map,
521 std::multimap<
unsigned int, std::pair<unsigned int, unsigned int>> &vertex_map )
523 if(vertex_crossings_map.size() == 0)
return;
526 double vdrift = 8.00;
527 double z_bunch_separation = 0.107 * vdrift;
529 for(
auto [ivert, crossing] : vertex_crossings_map)
536 if(
Verbosity() > 1) std::cout <<
" ivert " << ivert <<
" crossing " << crossing << std::endl;
538 if(crossing == 0)
continue;
542 auto ret = vertex_map.equal_range(ivert);
543 for(
auto it = ret.first;
it != ret.second; ++
it)
545 unsigned int tpcid =
it->second.first;
547 if(
Verbosity() > 1) std::cout <<
" tpcid " << tpcid <<
" original z " << tpc_track->
get_z() << std::endl;
565 if(
Verbosity() > 2) std::cout <<
" original local cluster z " << tpc_clus->
getLocalY() << std::endl;
569 double clus_z = global[2];
570 if(
Verbosity() > 2) std::cout <<
" global: x " << global[0] <<
" y " << global[1] <<
" z " << global[2] << std::endl;
572 double corrected_clus_z;
574 corrected_clus_z = clus_z + crossing * z_bunch_separation;
576 corrected_clus_z = clus_z - crossing * z_bunch_separation;
580 double surfZCenter = center[2];
581 double local_z = corrected_clus_z - surfZCenter;
583 if(
Verbosity() > 2) std::cout <<
" original clus_z " << clus_z <<
" corrected_clus_z " << corrected_clus_z <<
" local z " << local_z << std::endl;
597 std::cout <<
" global check: x " << global_check[0] <<
" y " << global_check[1] <<
" z " << global_check[2] << std::endl;
607 std::map<unsigned int, double> &vertex_crossings_map,
608 std::multimap<
unsigned int, std::pair<unsigned int, unsigned int>> &vertex_map,
609 std::map<unsigned int, int> &tpc_crossing_map )
611 if(vertex_crossings_map.size() == 0)
return;
615 std::multimap<unsigned int, std::pair<unsigned int, unsigned int>> bad_map;
616 for(
auto [ivert, crossing] : vertex_crossings_map)
618 if(
Verbosity() > 1) std::cout <<
" CleanVertexMap: ivert " << ivert <<
" crossing " << crossing << std::endl;
622 auto ret = vertex_map.equal_range(ivert);
623 for(
auto it = ret.first;
it != ret.second; ++
it)
625 unsigned int tpcid =
it->second.first;
626 unsigned int si_id =
it->second.second;
629 auto sit = tpc_crossing_map.find(tpcid);
630 if(
abs(sit->second - crossing) > 2)
634 std::cout <<
" Crossing mismatch: ivert " << ivert <<
" tpc id " << tpcid <<
" vert crossing " << crossing <<
" track crossing " << sit->second << std::endl;
636 bad_map.insert(std::make_pair(ivert, std::make_pair(tpcid, si_id)));
648 for(
auto [ivert, id_pair] : bad_map)
650 unsigned int tpc_id = id_pair.first;
651 unsigned int si_id = id_pair.second;
656 auto ret = vertex_map.equal_range(ivert);
657 for(
auto it = ret.first;
it != ret.second; ++
it)
659 if(
it->second.first == tpc_id &&
it->second.second == si_id)
661 if(
Verbosity() > 1) std::cout <<
" erasing entry for ivert " << ivert <<
" tpc_id " << tpc_id <<
" si_id " << si_id << std::endl;
662 vertex_map.erase(
it);
672 std::vector<double> &vertex_list,
673 std::multimap<
unsigned int, std::pair<unsigned int, unsigned int>> &vertex_map,
674 std::map<unsigned int, double> &vertex_crossings_map)
676 if(vertex_list.size() == 0)
return;
678 for(
unsigned int ivert = 0; ivert<vertex_list.size(); ++ivert)
680 std::vector<double> crossing_vec;
682 double vert_z = vertex_list[ivert];
683 if(
Verbosity() > 1) std::cout <<
"Vertex " << ivert <<
" vertex z " << vert_z << std::endl;
685 auto ret = vertex_map.equal_range(ivert);
686 for(
auto it = ret.first;
it != ret.second; ++
it)
688 unsigned int tpcid =
it->second.first;
690 double z_mismatch = tpc_z - vert_z;
693 crossing_vec.push_back(crossings);
697 <<
" tpc z " << tpc_z <<
" crossings " << crossings << std::endl;
700 if(crossing_vec.size() == 0)
continue;
701 double crossing_median =
getMedian(crossing_vec);
703 double crossing_avge = 0.0;
704 double crossing_wt = 0.0;
705 for(
auto cross : crossing_vec)
707 if(fabs(
cross-crossing_median) < 1.0)
709 crossing_avge +=
cross;
713 crossing_avge /= crossing_wt;
714 double crossing_avge_rounded =
round(crossing_avge);
715 vertex_crossings_map.insert(std::make_pair(ivert, crossing_avge_rounded));
717 if(
Verbosity() > 1) std::cout <<
" crossing_median " << crossing_median <<
" crossing average = " << crossing_avge <<
" crossing_wt " << crossing_wt
718 <<
" crossing integer " << crossing_avge_rounded << std::endl;
724 std::multimap<
double, std::pair<unsigned int, unsigned int>> &si_sorted_map,
725 std::vector<double> &vertex_list,
726 std::multimap<
unsigned int, std::pair<unsigned int, unsigned int>> &vertex_map)
728 if(si_sorted_map.size() == 0)
return;
734 unsigned int nvert = 0;
735 for(
auto it = si_sorted_map.begin();
it != si_sorted_map.end(); ++
it)
737 double z_si =
it->first;
738 std::pair<unsigned int, unsigned int> id_pair =
it->second;
739 if(
Verbosity() > 1) std::cout <<
" z_si " << z_si <<
" tpc_id " <<
it->second.first <<
" si_id " <<
it->second.second << std::endl;
743 vertex_map.insert(std::make_pair(nvert, id_pair));
751 vertex_list.push_back(zavge);
758 vertex_map.insert(std::make_pair(nvert, id_pair));
763 vertex_list.push_back(zavge);
769 std::set<unsigned int> &tpc_matched_set,
770 std::multimap<unsigned int, unsigned int> &tpc_matches )
785 <<
": Processing seed itrack: " << phtrk_iter->first
807 if(
_field.find(
"2d") != std::string::npos)
809 sign_phi_correction *= -1;
811 sign_phi_correction *= -1;
818 if(tpc_pt < 6.0) mag = 2;
819 if(tpc_pt < 3.0) mag = 4.0;
823 cout <<
"Original TPC tracklet:" << endl;
833 tpc_phi -=
fscdphi->Eval(tpc_eta);
836 if(tpc_phi < -
M_PI) tpc_phi += 2.0*
M_PI;
837 if(tpc_phi >
M_PI) tpc_phi -= 2.0*
M_PI;
861 cout <<
" tpc_phi " << tpc_phi <<
" si_phi " << si_phi <<
" dphi " << tpc_phi-si_phi <<
" phi search " <<
_phi_search_win*mag <<
" tpc_eta " << tpc_eta
862 <<
" si_eta " << si_eta <<
" deta " << tpc_eta-si_eta <<
" eta search " <<
_eta_search_win*mag << endl;
863 std::cout <<
" tpc x " << tpc_x <<
" si x " << si_x <<
" tpc y " << tpc_y <<
" si y " << si_y <<
" tpc_z " << tpc_z <<
" si z " << si_z << std::endl;
867 bool eta_match =
false;
868 bool phi_match =
false;
869 bool position_match =
false;
885 cout <<
" phi_search_win_lo " << phi_search_win_lo <<
" phi_search_win_hi " << phi_search_win_hi << endl;
887 if( (tpc_phi - si_phi) > phi_search_win_lo && (tpc_phi - si_phi) < phi_search_win_hi) phi_match =
true;
896 position_match =
true;
905 position_match =
true;
908 if(eta_match && phi_match && position_match)
918 cout <<
" tpc_phi " << tpc_phi <<
" si_phi " << si_phi <<
" phi_match " << phi_match
919 <<
" tpc_eta " << tpc_eta <<
" si_eta " << si_eta <<
" eta_match " << eta_match << endl;
920 std::cout <<
" tpc x " << tpc_x <<
" si x " << si_x <<
" tpc y " << tpc_y <<
" si y " << si_y <<
" tpc_z " << tpc_z <<
" si z " << si_z << std::endl;
925 cout <<
" Try_silicon: pt " << tpc_pt <<
" tpc_phi " << tpc_phi <<
" si_phi " << si_phi <<
" dphi " << tpc_phi-si_phi
926 <<
" tpc_eta " << tpc_eta <<
" si_eta " << si_eta <<
" deta " << tpc_eta-si_eta <<
" tpc_x " << tpc_x <<
" tpc_y " << tpc_y <<
" tpc_z " << tpc_z
927 <<
" dx " << tpc_x - si_x <<
" dy " << tpc_y - si_y <<
" dz " << tpc_z - si_z
935 std::multimap<unsigned int, unsigned int> additional_tpc_matches;
936 std::multimap<unsigned int, unsigned int> remove_tpc_matches;
937 std::set<unsigned int> additional_tpc_matched_set;
938 for(
unsigned int tpcid : tpc_matched_set)
940 auto ret = tpc_matches.equal_range(tpcid);
942 unsigned int size = std::distance(ret.first, ret.second);
943 if(
Verbosity() > 1) std::cout <<
" tpcid " << tpcid <<
" number of matches " << size << std::endl;
945 if(size == 1)
continue;
949 for(
auto it = ret.first;
it != ret.second; ++
it)
952 if(counter == 0)
continue;
954 unsigned int si_id =
it->second;
956 auto newTrack = std::make_unique<SvtxTrack_v2>();
958 if(
Verbosity() > 1) cout <<
"Extra match, add a new track to node tree with key " << lastTrackKey + 1 << endl;
960 newTrack->set_id(lastTrackKey+1);
962 newTrack->set_x(this_track->get_x());
963 newTrack->set_y(this_track->get_y());
964 newTrack->set_z(this_track->get_z());
966 newTrack->set_charge(this_track->get_charge());
967 newTrack->set_px(this_track->get_px());
968 newTrack->set_py(this_track->get_py());
969 newTrack->set_pz(this_track->get_pz());
971 for(
int i = 0; i < 6; ++i)
973 for(
int j = 0; j < 6; ++j)
975 newTrack->set_error(i,j, this_track->get_error(i,j));
983 iter != this_track->end_cluster_keys();
990 newTrack->insert_cluster_key(cluster_key);
995 if(
Verbosity() > 1) cout <<
" -- inserting new track with id " << newTrack->get_id() <<
" from TPC tracklet " << tpcid <<
" into trackmap " << endl;
999 remove_tpc_matches.insert(std::make_pair(tpcid, si_id));
1001 additional_tpc_matches.insert(std::make_pair(newTrack->get_id(), si_id));
1002 additional_tpc_matched_set.insert(newTrack->get_id());
1005 for(
auto [tpcid, si_id] : additional_tpc_matches)
1007 tpc_matched_set.insert(tpcid);
1008 tpc_matches.insert(std::make_pair(tpcid, si_id));
1009 if(
Verbosity() > 1) std::cout <<
"add to tpc_matches tpc_id " << tpcid <<
" si_id " << si_id << std::endl;
1011 for(
auto [tpcid, si_id] : remove_tpc_matches)
1013 if(
Verbosity() > 1) std::cout <<
"remove from tpc_matches tpc_id " << tpcid <<
" si_id " << si_id << std::endl;
1016 auto ret = tpc_matches.equal_range(tpcid);
1017 for(
auto it = ret.first;
it != ret.second; ++
it)
1019 if(
it->first == tpcid &&
it->second == si_id)
1021 tpc_matches.erase(
it);
1031 std::multimap<unsigned int, unsigned int> &tpc_matches,
1032 std::set<int> &crossing_set,
1033 std::multimap<
int, std::pair<unsigned int, unsigned int>> &crossing_matches,
1034 std::map<unsigned int, int> &tpc_crossing_map )
1037 for(
auto [tpcid, si_id] : tpc_matches)
1040 double tpc_z = tpc_track->
get_z();
1041 double tpc_pt = tpc_track->
get_pt();
1044 double si_z = si_track->
get_z();
1050 if(tpc_pt < 6.0) mag = 2;
1051 if(tpc_pt < 3.0) mag = 4.0;
1058 crossing_matches.insert(std::make_pair(crossing,std::make_pair(tpc_track->
get_id(), si_track->
get_id())));
1059 crossing_set.insert(crossing);
1060 tpc_crossing_map.insert(std::make_pair(tpc_track->
get_id(), crossing));
1062 std::cout <<
" triggered: tpc_trackid " << tpc_track->
get_id()
1063 <<
" eta " << tpc_track->
get_eta()
1064 <<
" pT " << tpc_track->
get_pt()
1065 <<
" si_trackid " << si_track->
get_id()
1066 <<
" z_tpc " << tpc_track->
get_z()
1067 <<
" z_si " << si_track->
get_z()
1068 <<
" crossing " << crossing
1076 std::multimap<unsigned int, unsigned int> &tpc_matches,
1077 std::set<int> &crossing_set,
1078 std::multimap<
int, std::pair<unsigned int, unsigned int>> &crossing_matches,
1079 std::map<unsigned int, int> &tpc_crossing_map )
1082 for(
auto [tpcid, si_id] : tpc_matches)
1085 double tpc_z = tpc_track->
get_z();
1089 double si_z = si_track->
get_z();
1093 crossing_matches.insert(std::make_pair(crossing,std::make_pair(tpc_track->
get_id(), si_track->
get_id())));
1094 crossing_set.insert(crossing);
1095 tpc_crossing_map.insert(std::make_pair(tpc_track->
get_id(), crossing));
1097 std::cout <<
" pileup: tpc_trackid " << tpc_track->
get_id()
1098 <<
" eta " << tpc_track->
get_eta()
1099 <<
" pT " << tpc_track->
get_pt()
1100 <<
" si_trackid " << si_track->
get_id()
1101 <<
" z_tpc " << tpc_track->
get_z()
1102 <<
" z_si " << si_track->
get_z()
1103 <<
" crossing " << crossing
1129 if( !cluster )
continue;