11 #include <HelixHough/HelixKalmanState.h>
12 #include <HelixHough/SimpleHit3D.h>
13 #include <HelixHough/SimpleTrack3D.h>
59 #define LogDebug(exp) std::cout<<"DEBUG: " <<__FILE__<<": "<<__LINE__<<": "<< exp
60 #define LogError(exp) std::cout<<"ERROR: " <<__FILE__<<": "<<__LINE__<<": "<< exp
61 #define LogWarning(exp) std::cout<<"WARNING: "<<__FILE__<<": "<<__LINE__<<": "<< exp
69 using namespace Eigen;
73 unsigned int min_nlayers,
76 _t_output_io(nullptr),
79 _min_nlayers(min_nlayers),
86 _use_max_kappa(
false),
87 fill_multi_zvtx(
false),
89 _max_kappa(numeric_limits<float>::
max()),
100 _ca_phi_cut(
M_PI/36.),
107 _override_min_zvtx_tracks(0),
128 _bbc_vertexes(nullptr),
131 _hough_space(nullptr),
132 _hough_funcs(nullptr),
133 _ntp_zvtx_by_event(nullptr),
134 _ntp_zvtx_by_track(nullptr),
138 _kappa_d_phi(nullptr),
141 _fname(
"init_zvtx_diagnostic_file.root"),
144 _layer_ilayer_map_all(),
159 separate_helicity(
false),
201 _ofile =
new TFile(
"z0_dzdl_kappa_phi_d.root",
"recreate");
211 for (
unsigned int izoom =0; izoom<
nzooms ; ++izoom)
245 _kappa_d_phi =
new TH3D(
"kappa_d_phi",
"kappa_d_phi",n_kappa_bins,
_hough_space->
get_kappa_min(),
_hough_space->
get_kappa_max(),n_d_bins,
_hough_space->
get_d_min(),
_hough_space->
get_d_max(),n_phi_bins,
_hough_space->
get_phi_min(),
_hough_space->
get_phi_max());
253 _ntp_zvtx_by_event =
new TNtuple(
"ntp_zvtx_by_event",
"all vertices found event-by-event",
"event:zvtx");
254 _ntp_zvtx_by_track =
new TNtuple(
"ntp_zvtx_by_track",
"track-by-track (zvtx + z0) distribution",
"event:zvtx");
263 <<
"====================== Leaving PHInitZVertexing::Setup() ======================"
265 cout <<
" Magnetic field set to: " <<
_mag_field <<
" Tesla" << endl;
266 cout <<
" Number of tracking layers: " <<
_nlayers << endl;
267 for (
unsigned int i = 0; i <
_nlayers; ++i) {
268 cout <<
" Tracking layer #" << i <<
" " <<
"radius = "
273 cout <<
" Minimum pT: " <<
_min_pt << endl;
274 cout <<
" Maximum DCA: " << boolalpha <<
_cut_on_dca << noboolalpha
277 cout <<
" Maximum DCA cut: " <<
_dcaxy_cut << endl;
279 cout <<
" Maximum DCAZ cut: " <<
_dcaz_cut << endl;
280 cout <<
" Momentum rescale factor: " <<
_pt_rescale << endl;
282 <<
"==========================================================================="
292 if (
Verbosity() > 1) cout <<
"PHInitZVertexing::Process -- entered" << endl;
327 bool zvtx_found =
false;
334 unsigned int nseq = 1;
335 unsigned int nattempt = 1;
341 if (iseq==nseq)
break;
354 for (
unsigned int iattempt =0; iattempt<nattempt; ++iattempt ){
366 unsigned int zoomlevel = 0;
371 cout<<
"InitZVertexing:: nkappa " <<
nkappa<<
" nphi "<<
nphi<<
" nd "<<
nd<<
" ndzdl "<<
ndzdl<<
" nz0 " <<
nz0<<endl;
378 while(it_begin != it_end)
380 delete it_begin->second;
385 while(it_begin != it_end)
387 delete it_begin->second;
392 while(it_begin != it_end)
394 delete it_begin->second;
399 while(it_begin != it_end)
401 delete it_begin->second;
421 for (zoomlevel =1; zoomlevel<2; ++zoomlevel){
439 cout<<
"CellularAutomaton failed. "<<endl;
453 if (code==-1) zvtx_found =
false;
463 }
else if (iseq>0 && code<0 && !zvtx_found){
470 cout<<
"z-vertex not fitted. "<<endl;
472 cout<<
"z-vertex not found. "<<endl;
475 if (iseq<nseq) cout<<
"z-vertex fitted. "<<endl;
526 std::vector<unsigned int> zoom {n_kappa, n_phi, n_d, n_dzdl, n_z0};
537 "PHCompositeNode",
"DST"));
539 cerr <<
PHWHERE <<
"DST Node missing, doing nothing." << endl;
549 "SvtxTrackMap_zvtx",
"PHObject");
551 if (
Verbosity() > 0) cout <<
"Svtx/SvtxTrackMap node added" << endl;
573 _radii.assign(_nlayers, 0.0);
574 map<float, int> radius_layer_map;
582 layerrange.first; layeriter != layerrange.second; ++layeriter) {
584 radius_layer_map.insert(
585 make_pair(layeriter->second->get_radius(),
586 layeriter->second->get_layer()));
592 laddergeos->get_begin_end();
594 layerrange.first; layeriter != layerrange.second; ++layeriter) {
596 radius_layer_map.insert(
597 make_pair(layeriter->second->get_radius(),
598 layeriter->second->get_layer()));
602 if (mvtxladdergeos) {
604 mvtxladdergeos->get_begin_end();
606 layerrange.first; layeriter != layerrange.second; ++layeriter) {
608 radius_layer_map.insert(
609 make_pair(layeriter->second->get_radius(),
610 layeriter->second->get_layer()));
619 for (map<float, int>::iterator iter = radius_layer_map.begin();
620 iter != radius_layer_map.end(); ++iter) {
628 if(
Verbosity() > 1) cout <<
" PHInitZVertexing: inserting _layer_ilayer_map entry number " <<ilayer <<
" for layer " << iter->second << endl;
640 for (; miter != begin_end.second; ++miter) {
662 laddergeos->get_begin_end();
664 for (; miter != begin_end.second; ++miter) {
682 if (mvtxladdergeos) {
684 mvtxladdergeos->get_begin_end();
686 for (; miter != begin_end.second; ++miter) {
708 map<int, float>::iterator mat_it;
724 unsigned int clusid = 0;
725 unsigned int ilayer = 0;
728 cout <<
PHWHERE <<
"Get clusters:" << endl;
729 cout <<
" _cluster_map has " <<
_cluster_map->
size() <<
" entries " << endl;
730 cout <<
" _layer_ilayer_map has " <<
_layer_ilayer_map.size() <<
" entries" << endl;
734 for (
auto hitsetitr = hitsetrange.first;
735 hitsetitr != hitsetrange.second;
738 for(
auto clusIter = range.first; clusIter != range.second; ++clusIter ){
752 if(
Verbosity() > 200) cout <<
" clusid " << clusid <<
" cluskey " << cluskey <<
" layer " << layer <<
" ilayer " << ilayer << endl;
765 for (
int i = 0; i < 3; ++i) {
766 for (
int j = i; j < 3; ++j) {
773 hits_map.insert(std::pair<unsigned int, SimpleHit3D>(hit3d.
get_id(),hit3d));
785 if(rough_tracknum < 200)
787 else if(rough_tracknum < 1000)
789 else if(rough_tracknum < 3000)
796 <<
" for this event, which has " <<
_nclus_mvtx <<
" MVTX clusters " << endl;
801 cout <<
"-------------------------------------------------------------------"
803 cout <<
"PHInitZVertexing::process_event has the following clusters copied to SimpleHit3D:"
806 cout <<
"n init clusters = " <<
hits_map.size() << endl;
808 for (std::map<unsigned int,SimpleHit3D>::iterator
it=
hits_map.begin();
817 cout <<
"-------------------------------------------------------------------"
832 std::vector<SvtxVertex_v1> svtx_vertex_list;
834 for (
unsigned int vid = 0; vid < nvertex; ++vid ){
839 for (
int i = 0; i < 3; ++i)
853 svtx_vertex_list.push_back(vertex);
860 std::vector<SimpleHit3D> track_hits;
864 for (
unsigned int itrack = 0; itrack <
_tracks.size(); ++itrack) {
865 if (
_tracks.at(itrack).vertex_id >= nvertex)
continue;
869 track_hits =
_tracks.at(itrack).hits;
871 for (
unsigned int ihit = 0;
ihit < track_hits.size(); ++
ihit) {
875 clusterkey = track_hits.at(
ihit).get_cluskey();
880 <<
": itrack: " << itrack
881 <<
": nhits: " << track_hits.size()
882 <<
": clusterID: " << clusterID
883 <<
": clusterkey: " << clusterkey
884 <<
": layer: " <<
layer
891 float kappa =
_tracks.at(itrack).kappa;
899 float x_center = cos(phi) * (d + 1 / kappa);
900 float y_center = sin(phi) * (d + 1 / kappa);
904 if ((track_hits[0].get_x()- x_center)*(track_hits[track_hits.size()-1].get_y()- y_center)
905 - (track_hits[0].get_y()- y_center)*(track_hits[track_hits.size()-1].get_x()- x_center)> 0)
914 pZ = pT * dzdl / sqrt(1.0 - dzdl * dzdl);
916 int ndf = 2 *
_tracks.at(itrack).hits.size() - 5;
919 track.
set_px(pT * cos(phi - helicity *
M_PI / 2));
920 track.
set_py(pT * sin(phi - helicity *
M_PI / 2));
939 for (
unsigned int row = 0; row < 6; ++row) {
940 for (
unsigned int col = 0;
col < 6; ++
col) {
945 unsigned int vid =
_tracks.at(itrack).vertex_id;
947 track.
set_x(svtx_vertex_list[vid].get_x() + d * cos(phi));
948 track.
set_y(svtx_vertex_list[vid].get_y() + d * sin(phi));
949 track.
set_z(svtx_vertex_list[vid].get_z() + z0);
952 svtx_vertex_list[vid].insert_track(track.
get_id());
955 cout <<
"track " << itrack <<
" quality = " << track.
get_quality()
957 cout <<
"px = " << track.
get_px() <<
" py = " << track.
get_py()
958 <<
" pz = " << track.
get_pz() << endl;
968 for (
unsigned int vid = 0; vid <
_vertex_list.size(); ++vid ){
1013 for (
id=0;
id<
nd; ++
id) {
1026 bins_map.insert(std::pair<unsigned int, HelixHoughBin*>(
bin,_bin));
1040 unsigned int cluster_id=-999;
1041 std::vector<float> kappa_phi_d_ranges;
1048 kappa_phi_d_ranges.push_back(kappa_min);
1049 kappa_phi_d_ranges.push_back(kappa_max);
1050 kappa_phi_d_ranges.push_back(phi_min);
1051 kappa_phi_d_ranges.push_back(phi_max);
1052 kappa_phi_d_ranges.push_back(d_min);
1053 kappa_phi_d_ranges.push_back(d_max);
1060 for (std::map<unsigned int,SimpleHit3D>::iterator
it=
hits_map.begin();
1065 cluster_id =
it->first;
1069 auto hitused =
hits_used.find(cluster_id);
1071 used =
hits_used.find(cluster_id)->second;
1074 hitpos3d[0] = hit.
get_x();
1075 hitpos3d[1] = hit.
get_y();
1076 hitpos3d[2] = hit.
get_z();
1081 std::vector<float> z0_range;
1084 z0_range.push_back(z0_min);
1085 z0_range.push_back(z0_max);
1086 float dzdl_range[2] = {4,5};
1096 unsigned int dzdl_bin_range[2];
1099 if (dzdl_range[0] > dzdl_max)
continue;
1100 else if (dzdl_range[0] <dzdl_min) dzdl_bin_range[0]=0;
1101 if (dzdl_range[1] < dzdl_min)
continue;
1104 cout<<
"dzdl_min " <<dzdl_range[0]<<
" dzdl_max "<<dzdl_range[1]<<endl;
1106 for (
il =dzdl_bin_range[0];
il<dzdl_bin_range[1]+1; ++
il) {
1107 unsigned int zbins[5] = {0,0,0,
il,
iz};
1110 cout<<
"bin "<<
bin<<endl;
1115 bins_map.find(
bin)->second->add_cluster_ID(cluster_id);
1118 cout<<
"count "<<
bins_map.find(
bin)->second->get_count()<<endl;
1124 kappa_phi_d_ranges.clear();
1125 hitpos3d[0]=-999; hitpos3d[1]= -999; hitpos3d[2]=-999;
1128 if(
Verbosity() > 1) cout<<
"total number of clusters "<<icluster<<endl;
1138 unsigned int max_i=-1;
1139 unsigned int max_count = 0;
1143 unsigned int bins[5]={0,0,0,j,i};
1145 unsigned int count =
bins_map.find(
bin)->second->get_count();
1146 if (count <1)
continue;
1147 if (count>max_count) {
1152 cout<<
"bin "<<
bin<<endl;
1153 cout<<
"z0 bin "<<i<<
" dzdl bin "<<j<<
" count "<< count<<endl;
1154 cout<<
"z0 bin "<<
bins_map.find(
bin)->second->get_z0_bin(0)<<
" dzdl bin "<<
bins_map.find(
bin)->second->get_dzdl_bin(0)<<endl;
1175 unsigned int bins[5]={0,0,0,
k,max_i};
1177 unsigned int count =
bins_map.find(
bin)->second->get_count();
1178 if (count <1)
continue;
1186 if(
Verbosity() > 1) cout<<
"bins_map_sel.size at zoom "<<zoomlevel <<
" : (find_track_candidates_z_init) " <<
bins_map_sel.size()<<endl;
1196 unsigned int cluster_id=-999;
1203 unsigned int izprev = houghbin->
get_z0_bin(zoomlevel-1);
1204 unsigned int ilprev = houghbin->
get_dzdl_bin(zoomlevel-1);
1205 unsigned int ipprev = houghbin->
get_phi_bin(zoomlevel-1);
1218 hitpos3d[0] = hit.
get_x();
1219 hitpos3d[1] = hit.
get_y();
1220 hitpos3d[2] = hit.
get_z();
1225 std::vector<float> kappa_phi_d_ranges;
1232 kappa_phi_d_ranges.push_back(kappa_min);
1233 kappa_phi_d_ranges.push_back(kappa_max);
1234 kappa_phi_d_ranges.push_back(phi_min);
1235 kappa_phi_d_ranges.push_back(phi_max);
1236 kappa_phi_d_ranges.push_back(d_min);
1237 kappa_phi_d_ranges.push_back(d_max);
1248 std::vector<float> z0_range;
1251 z0_range.push_back(z0_min);
1252 z0_range.push_back(z0_max);
1253 float dzdl_range[2];
1257 unsigned int dzdl_bin_range[2];
1264 if (dzdl_range[0] > dzdl_max)
continue;
1265 else if (dzdl_range[0] <dzdl_min) dzdl_bin_range[0]=0;
1266 if (dzdl_range[1] < dzdl_min)
continue;
1271 for (
il =dzdl_bin_range[0];
il<dzdl_bin_range[1]+1; ++
il) {
1272 unsigned int zbins[5] = {0,0,0,
il,
iz};
1274 houghbin->
set_bin(zoomlevel,curbin);
1275 houghbin->
set_bins(zoomlevel,curbin);
1283 bins_map_cur.find(curgbin)->second->add_cluster_ID(cluster_id);
1291 cur_bin->
set_bin(zoomlevel,curbin);
1292 cur_bin->
set_bins(zoomlevel,curbin);
1294 bins_map_cur.insert(std::pair<unsigned int,HelixHoughBin*>(curgbin,cur_bin));
1303 unsigned int count =
bins_map_cur.find(curgbin)->second->get_count();
1311 hitpos3d[0]=-999; hitpos3d[1]= -999; hitpos3d[2]=-999;
1312 kappa_phi_d_ranges.clear();
1318 if(
Verbosity() > 1) cout<<
"bins_map_cur.size at zoom "<<zoomlevel <<
" (vote_z) : " <<
bins_map_cur.size()<<endl;
1331 unsigned int count = houghbin->
get_count();
1332 if (count <1)
continue;
1349 cout<<
"il "<<
il<<
" iz "<<
iz<<
" count "<< count<<endl;
1357 if(
Verbosity() > 1) cout<<
"bins_map_sel.size at zoom "<< zoomlevel<<
" (find_track_candidates_z) : " <<
bins_map_sel.size()<<endl;
1364 unsigned int ilfill= -1;
1365 unsigned int izfill= -1;
1371 unsigned int phi_bin_range[2];
1381 if (count==0){ ilfill =
il; izfill =
iz;}
1383 unsigned int ikprev =0;
1384 unsigned int ipprev =0;
1385 unsigned int idprev =0;
1389 idprev = houghbin->
get_d_bin(zoomlevel-1);
1392 cout<<
"bin "<<
bin<<endl;
1393 cout<<
"il " <<
il<<
" iz "<<
iz<<
" ikprev "<<ikprev<<
" ipprev "<<ipprev<<
" idprev "<<idprev<<
" count"<<houghbin->
get_count()<<endl;
1410 unsigned int cluster_id = *iter;
1416 cout<<
"cluster id "<<cluster_id<<
" x "<<
hitpos3d[0]<<
" y "<<
hitpos3d[1]<<endl;
1419 for (
id=0;
id<
nd; ++
id){
1420 std::vector<float> kappa_d_ranges;
1434 kappa_d_ranges.push_back(kappa_min);
1435 kappa_d_ranges.push_back(kappa_max);
1436 kappa_d_ranges.push_back(d_min);
1437 kappa_d_ranges.push_back(d_max);
1441 float phi_r_range[2];
1442 float phi_l_range[2];
1443 unsigned int phi_r_bin_range[2];
1444 unsigned int phi_l_bin_range[2];
1455 bool phi_r_out_of_range=
false;
1456 bool phi_l_out_of_range=
false;
1457 if (phi_r_range[0] > phi_max) phi_r_out_of_range=
true;
1458 else if (phi_r_range[0] <phi_min) phi_r_bin_range[0]=0;
1459 if (phi_r_range[1] < phi_min) phi_r_out_of_range=
true;
1462 if (phi_l_range[0] > phi_max) phi_r_out_of_range=
true;
1463 else if (phi_r_range[0] <phi_min) phi_r_bin_range[0]=0;
1464 if (phi_r_range[1] < phi_min) phi_r_out_of_range=
true;
1467 for (
int i = 0; i<2;++i){
1469 if ((i==0&&phi_r_out_of_range) || (i==1&&phi_l_out_of_range))
continue;
1472 phi_bin_range[0]=phi_r_bin_range[0];
1473 phi_bin_range[1]=phi_r_bin_range[1];
1476 phi_bin_range[0]=phi_l_bin_range[0];
1477 phi_bin_range[1]=phi_l_bin_range[1];
1482 for (
ip= phi_bin_range[0];
ip<phi_bin_range[1]+1;++
ip){
1485 houghbin->
set_bin(zoomlevel,curbin);
1486 houghbin->
set_bins(zoomlevel,curbin);
1492 bins_map_cur.find(curgbin)->second->add_cluster_ID(cluster_id);
1500 cur_bin->
set_bin(zoomlevel,curbin);
1501 cur_bin->
set_bins(zoomlevel,curbin);
1503 bins_map_cur.insert(std::pair<unsigned int,HelixHoughBin*>(curgbin,cur_bin));
1514 kappa_d_phi_map.insert(std::pair<unsigned int,unsigned int>(curgbin,1));
1517 unsigned int kpbins[5]={
ik,
ip,0,
il,
iz};
1526 kappa_phi_map.insert(std::pair<unsigned int,unsigned int>(kpbin,1));
1529 unsigned int dpbins[5]={0,
ip,
id,
il,
iz};
1538 d_phi_map.insert(std::pair<unsigned int,unsigned int>(dpbin,1));
1542 if (
il ==ilfill &&
iz == izfill){
1549 cout <<
"il "<<
il<<
" iz "<<
iz<<endl;
1562 float phi_prev_range[2];
1563 float phi_next_range[2];
1564 unsigned int phi_bin_range[2];
1572 phi_prev_range[0]= phi_next_range[0];
1573 phi_prev_range[1]= phi_next_range[1];
1580 bool phi_out_of_range=
false;
1581 if (phi_range[0] > phi_max) phi_out_of_range=
true;
1582 else if (phi_range[0] <phi_min) phi_bin_range[0]=0;
1583 if (phi_range[1] < phi_min) phi_out_of_range=
true;
1586 if (phi_out_of_range)
continue;
1588 for (
ip= phi_bin_range[0];
ip<phi_bin_range[1]+1;++
ip){
1591 houghbin->
set_bin(zoomlevel,curbin);
1592 houghbin->
set_bins(zoomlevel,curbin);
1595 cout<<
" curgbin "<<curgbin<<endl;
1600 bins_map_cur.find(curgbin)->second->add_cluster_ID(cluster_id);
1608 cur_bin->
set_bin(zoomlevel,curbin);
1609 cur_bin->
set_bins(zoomlevel,curbin);
1611 bins_map_cur.insert(std::pair<unsigned int,HelixHoughBin*>(curgbin,cur_bin));
1613 cout<<
" il "<<
il<<
" iz "<<
iz<<
" il from bin"<<
bins_map_cur.find(curbin)->second->get_dzdl_bin(0)<<
" iz from bin "<<
bins_map_cur.find(curbin)->second->get_z0_bin(0)<<endl;
1636 kappa_d_ranges.clear();
1642 if(
Verbosity() > 1) cout<<
"bins_map_cur.size at zoom "<<zoomlevel <<
" (vote_xy) : " <<
bins_map_cur.size()<<endl;
1657 unsigned int count = houghbin->
get_count();
1658 if (count <1)
continue;
1667 cout<<
"bin "<<
bin<<endl;
1668 cout<<
"kappa bin "<<
ik<<
" phi bin "<<
ip<<
" d bin "<<
id <<
" il "<<
il<<
" iz "<<
iz<<
" count "<< count<<endl;
1676 cout<<
"inserting selected bin.. "<<endl;
1677 cout<<
"bin "<<
bin<<endl;
1678 cout<<
"kappa bin "<<
ik<<
" phi bin "<<
ip<<
" d bin "<<
id <<
" il "<<
il<<
" iz "<<
iz<<
" count "<< count<<endl;
1685 if(
Verbosity() > 1) cout<<
"bins_map_prev.size at zoom "<<zoomlevel <<
" (find_track_candidates_xy) " <<
bins_map_prev.size()<<endl;
1695 unsigned int nclust = houghbin->
get_count();
1704 unsigned int var[3] = {1<<3,1<<4,(1<<3)+(1<<4)};
1705 unsigned int sign[3][4] = {{0<<3,1<<3,0,0},
1707 {(0<<3)+(0<<4),(0<<3)+(1<<4),(1<<3)+(0<<4),(1<<3)+(1<<4)}};
1708 for (
unsigned int lz =0; lz<3; ++lz){
1709 for (
unsigned int ss = 0; ss<4; ++ss ){
1710 if (lz<2 && ss>1)
continue;
1713 cout<<
"var "<<var[lz]<<
" ss "<<ss<<
" neighbor_global bin "<<neighborbin<<endl;
1715 if (globalbin==neighborbin)
continue;
1719 cout<<
"neighbor bin found"<<endl;
1721 unsigned int nclust_neighbor =
bins_map_sel.find(neighborbin)->second->get_count();
1722 if(nclust >= nclust_neighbor && nclust <= 1.1*
_nlayers) {
1724 cout<<
"merging "<<endl;
1727 iter !=
bins_map_sel.find(neighborbin)->second->end_clusters();
1741 if(
Verbosity() > 1) cout<<
"bins_map_sel.size at zoom " <<zoomlevel<<
" (prune_z) : " <<
bins_map_sel.size()<<endl;
1747 cout<<
"bins_map_prev.size at zoom " <<zoomlevel<<
" (pre prune_xy) : " <<
bins_map_prev.size()<<endl;
1750 unsigned int nclust = houghbin->
get_count();
1753 cout<<
"global bin "<<globalbin<<endl;
1760 unsigned int var[6] = {1,1<<1,1<<2,(1)+(1<<1),(1)+(1<<2),(1<<1)+(1<<2)};
1761 unsigned int sign[6][4]={{0,1,0,0},
1764 {0,1<<1,1,(1)+(1<<1)},
1765 {0,1<<2,1, (1)+(1<<2)},
1766 {0,1<<2,1<<1,(1<<1)+(1<<2)}};
1767 for (
unsigned int kpd =0; kpd<6; ++kpd){
1768 for (
unsigned int ss = 0; ss<4; ++ss ){
1769 if (kpd<3 && ss>1)
continue;
1772 if (globalbin==neighborbin)
continue;
1776 unsigned int nclust_neighbor =
bins_map_prev.find(neighborbin)->second->get_count();
1777 if(nclust >= nclust_neighbor && nclust <= 1.1*
_nlayers) {
1779 cout<<
"merging "<<endl;
1782 iter !=
bins_map_prev.find(neighborbin)->second->end_clusters();
1796 if(
Verbosity() > 1) cout<<
"bins_map_prev.size at zoom " <<zoomlevel<<
" (prune_xy) : " <<
bins_map_prev.size()<<endl;
1802 for (std::map<unsigned int,bool>::iterator
it =
hits_used.begin();
1810 unsigned int icluster = 0;
1822 track.
hits.push_back(hit);
1823 track.
hits[icluster].set_id(*iter);
1831 cout<<
"bins_map_prev selected, total bins : "<<
bins_map_prev.size()<<endl;
1847 auto search =
hits_map.find(*iter);
1851 track.
hits[icluster] = hit;
1852 track.
hits[icluster].set_id(*iter);
1862 new_tracks.push_back(track);
1872 unsigned int layer0 = 0;
1873 unsigned int layer1 = 1;
1874 unsigned int layer2 = 2;
1881 std::vector<unsigned int> layer_sorted_0;
1882 std::vector<unsigned int> layer_sorted_1;
1883 std::vector<unsigned int> layer_sorted_2;
1885 for (std::map<unsigned int,SimpleHit3D>::iterator
it=
hits_map.begin();
1889 unsigned int cluster_id =
it->first;
1891 auto hitused =
hits_used.find(cluster_id);
1893 used =
hits_used.find(cluster_id)->second;
1898 unsigned int hitid = hit.
get_id();
1900 if(
Verbosity() > 0) cout<<
"layer "<< layer<<
" hitid "<<hitid<<endl;
1901 if (layer == layer0 ) layer_sorted_0.push_back(hitid);
1902 else if (layer == layer1) layer_sorted_1.push_back(hitid);
1903 else if (layer == layer2) layer_sorted_2.push_back(hitid);
1909 cout<<
"layer 0 " <<layer_sorted_0.size()<<endl;
1910 cout<<
"layer 1 " <<layer_sorted_1.size()<<endl;
1911 cout<<
"layer 2 " <<layer_sorted_2.size()<<endl;
1914 std::set<unsigned int> clusters;
1915 bool fill_track =
false;
1916 unsigned int cluster_layer0 = 0;
1917 for (std::vector<unsigned int>::iterator it0 = layer_sorted_0.begin() ; it0 != layer_sorted_0.end(); ++it0)
1920 cout<<
"icluster layer 0 "<<cluster_layer0<<endl;
1925 if(
Verbosity() > 1) cout<<
"cluster_id for hit 0 "<< *it0<<endl;
1926 auto search0 =
hits_map.find(*it0);
1927 if (search0 ==
hits_map.end())
continue;
1929 float phi_h0 = atan2(hit0.
get_y(),hit0.
get_x());
1931 float z_h0 = hit0.
get_z();
1935 clusters.insert(*it0);
1937 for (std::vector<unsigned int>::iterator it1 = layer_sorted_1.begin(); it1 != layer_sorted_1.end(); ++it1 )
1940 float phi_h1 = atan2(hit1.
get_y(),hit1.
get_x());
1942 float z_h1 = hit1.
get_z();
1944 float phi_h0h1 = phi_h1-phi_h0;
1945 if (phi_h1< M_PI/2. && phi_h0 > 3*
M_PI/2.) phi_h0h1 += 2.*
M_PI;
1946 else if (phi_h1> 3*
M_PI/2. && phi_h0 <
M_PI/2.) phi_h0h1 -= 2.*
M_PI;
1951 for(std::vector<unsigned int>::iterator it2 = layer_sorted_2.begin(); it2 != layer_sorted_2.end(); ++it2)
1954 float phi_h2 = atan2(hit2.
get_y(),hit2.
get_x());
1956 float phi_h1h2 = phi_h2-phi_h1;
1957 if (phi_h2< M_PI/2. && phi_h1 > 3*
M_PI/2.) phi_h1h2 += 2.*
M_PI;
1958 else if (phi_h2> 3*
M_PI/2. && phi_h1 <
M_PI/2.) phi_h1h2 -= 2.*
M_PI;
1961 float z_h2 = hit2.
get_z();
1962 if (fabs(phi_h1h2) <
_ca_phi_cut && (z_h1-z_h0)*(z_h2-z_h1)>0 && fabs(z_h2-z_h1)<
_z_cut ){
1964 clusters.insert(*it1);
1965 clusters.insert(*it2);
1972 if(
Verbosity() > 1) cout<<
" fill_track "<<endl;
1973 unsigned int nclusters =0 ;
1975 for (std::set<unsigned int>::iterator
it=clusters.begin();
it!=clusters.end(); ++
it ){
1976 if(
Verbosity() > 1) cout<<
"cluster id "<<*
it<<endl;
1978 track.
hits[nclusters] = hit;
1979 track.
hits[nclusters].set_id(*
it);
1982 new_tracks.push_back(track);
1987 if(
Verbosity() > 1) cout<<
"number of triplets "<<new_tracks.size()<<endl;
1994 if(
Verbosity() > 1) cout<<
"Entering cellular autumaton : processing "<< candidate_tracks.size()<<
" tracks. "<<endl;
2018 for(
unsigned int i=0; i<candidate_tracks.size(); ++i) candidate_tracks[i].reset();
2019 candidate_tracks.clear();
2080 if(
Verbosity() > 1) cout<<
"Enter fit_vertex: all tracks vector size " <<
_tracks.size()<<endl;
2082 if (
_tracks.empty())
return -1;
2083 std::vector<SimpleTrack3D> vtx_tracks;
2087 unsigned int nzvtx = vtx_tracks.size();
2090 unsigned int nzbins= 2*1280;
2092 float zvalues[2*1280];
2093 unsigned int zcounts[2*1280];
2094 for (
unsigned int i = 0; i<nzbins; ++i)
2096 zvalues[i] =
_min_z0 + (i+0.5)*binsize;
2101 for (
unsigned int i = 0; i < nzvtx; ++i)
2103 double zvtx = vtx_tracks[i].z0;
2104 if (zvtx != zvtx || zvtx<_min_z0 || zvtx>
_max_z0)
continue;
2105 unsigned int zbin = (zvtx-
_min_z0)/binsize;
2110 std::vector<float> zvertices;
2111 std::vector<float> zmeans;
2112 std::vector<float> zsums;
2113 std::vector<float> zsigmas;
2114 std::vector<unsigned int> nzvtxes;
2116 for (
unsigned int j=2; j<(nzbins-3); ++j)
2118 if(
Verbosity() > 200 && zcounts[j] > 0) cout<<
"z countz "<<j<<
" "<<zcounts[j]<<endl;
2123 bool twobinspeak =
_mult_twobins*(zcounts[j-1]+zcounts[j-2])<(zcounts[j]+zcounts[j+1])
2124 &&
_mult_twobins*(zcounts[j+2]+zcounts[j+3])< (zcounts[j]+zcounts[j+1]);
2130 float zvertex=-999.;
2135 if (onebinpeak) zvertex = zvalues[j];
2136 if (twobinspeak) zvertex = (zvalues[j]*zcounts[j] + zvalues[j+1]*zcounts[j+1])/(zcounts[j]+zcounts[j+1]);
2137 zvertices.push_back(zvertex);
2138 nzvtxes.push_back(0);
2139 zmeans.push_back(0.0);
2140 zsums.push_back(0.0);
2141 zsigmas.push_back(0.0);
2143 if(
Verbosity() > 1) cout<<
" zbin " << j <<
" added one bin vertex at " << zvertex <<
" to nvertices with " << zcounts[j] <<
" contributing tracks " << endl;
2145 if(
Verbosity() > 1) cout<<
" zbin " << j <<
" added two bin vertex at " << zvertex <<
" to nvertices with " << zcounts[j]+zcounts[j+1] <<
" contributing tracks " << endl;
2151 cout<<
"number of candidate z-vertices : "<< zvertices.size()<<endl;
2152 for (
unsigned int j = 0; j <zvertices.size(); ++j)
2153 cout<<
"vertex " << j <<
" primary vertex position : "<<zvertices[j]<<endl;
2156 unsigned int izvtx= 999;
2157 for (
unsigned int i = 0; i < nzvtx; ++i) {
2158 double zvtx = vtx_tracks[i].z0;
2159 if (zvtx != zvtx)
continue;
2160 bool pass_dcaz=
false;
2161 for (
unsigned int k = 0;
k<zvertices.size(); ++
k) {
2162 bool new_vtx = fabs(zvtx-zvertices[
k]) <
_dcaz_cut;
2163 if (new_vtx) izvtx=
k;
2164 pass_dcaz = pass_dcaz || new_vtx;
2167 if (!pass_dcaz || izvtx>99)
continue;
2169 zsums[izvtx] += zvtx;
2172 if(
Verbosity() > 1) cout <<
"Loop over vertices and check for tracks that pass the dcaz cut of " <<
_dcaz_cut << endl;
2173 for (
unsigned int iz = 0;
iz<zvertices.size();++
iz ){
2174 if(
Verbosity() > 1) cout <<
" vertex " <<
iz <<
" ntracks passed " << nzvtxes[
iz] << endl;
2175 if (nzvtxes[
iz]==0)
continue;
2176 zmeans[
iz] = zsums[
iz]/nzvtxes[
iz];
2177 if(
Verbosity() > 1) cout<<
"zmean for passed vertex "<<
iz <<
" = "<<zmeans[
iz]<<endl;
2178 zvertices[
iz] = zmeans[
iz];
2183 ntp_data[1] = zmeans[
iz];
2188 for (
unsigned int j = 0; j < nzvtx; ++j){
2189 double zvtx = vtx_tracks[j].z0;
2190 if (zvtx != zvtx)
continue;
2196 ntp_data[1] = vtx_tracks[j].z0;
2202 bool pass_dcaz=
false;
2203 for (
unsigned int k = 0;
k<zvertices.size(); ++
k) {
2204 bool new_vtx = fabs(zvtx-zvertices[
k]) <
_dcaz_cut;
2205 if (new_vtx) izvtx=
k;
2206 pass_dcaz = pass_dcaz || new_vtx;
2209 if (!pass_dcaz || izvtx>99)
continue;
2210 zsums[izvtx] += pow(fabs(zmeans[izvtx] - vtx_tracks[j].
z0),2);
2213 std::vector<float> _multi_vtx;
2214 std::vector<std::vector<SimpleTrack3D> > _multi_vtx_tracks;
2215 std::vector<std::vector<double> > _multi_vtx_track_errors;
2216 std::vector<std::vector<Eigen::Matrix<float, 5, 5> > > _multi_vtx_track_covars;
2219 for (
unsigned int i = 0; i < zvertices.size(); ++i)
2221 if (nzvtxes[i]==0)
continue;
2222 _multi_vtx.push_back(zvertices[i]);
2223 std::vector<SimpleTrack3D> one_vtx_tracks;
2224 _multi_vtx_tracks.push_back(one_vtx_tracks);
2225 std::vector<double> one_track_errors;
2226 _multi_vtx_track_errors.push_back(one_track_errors);
2227 std::vector<Eigen::Matrix<float, 5, 5> > one_track_covars;
2228 _multi_vtx_track_covars.push_back(one_track_covars);
2229 zsigmas[i] = sqrt(zsums[i]/(nzvtxes[i]-1));
2230 if(
Verbosity() > 1) cout<<
" Add passed vertex " << i <<
" to _multi_vtx, with zvertex " << zvertices[i] <<
" and zsigma = "<<zsigmas[i]<<endl;
2233 unsigned int nzvtx_final = 0;
2234 for (
unsigned int k = 0;
k < nzvtx; ++
k){
2236 float zvtx = vtx_tracks[
k].z0;
2237 bool pass_dcaz=
false;
2238 for (
unsigned int iz = 0;
iz<_multi_vtx.size(); ++
iz) {
2239 bool new_vtx = fabs(zvtx-_multi_vtx[
iz]) <
_dcaz_cut;
2240 if (new_vtx) izvtx=
iz;
2241 pass_dcaz = pass_dcaz || new_vtx;
2244 if (!pass_dcaz || izvtx >= _multi_vtx.size())
continue;
2245 if(
Verbosity() > 1) cout<<
"adding a track to _multi_vtx number "<<izvtx<<endl;
2249 _multi_vtx_tracks[izvtx].push_back(vtx_tracks[
k]);
2251 _multi_vtx_track_errors[izvtx].push_back(
_track_states[k].chi2);
2255 if(
Verbosity() > 1) cout<<
"start final fitting of vertices in _multi_vtx.. "<<endl;
2256 for (
unsigned int i = 0; i<_multi_vtx.size(); ++i)
2258 if(
Verbosity() > 1) cout <<
" _multi_vtx " << i <<
" has " << _multi_vtx_tracks[i].size() <<
" tracks " << endl;
2259 if (_multi_vtx_tracks[i].size()==0)
continue;
2262 cout <<
" seed track vertex pre-fit: "
2278 cout<<
" - number of fitted tracks for vertex "<<i<<
" : "<<
n_vtx_tracks <<endl;
2279 cout <<
" seed track vertex post-fit: "
2292 cout<<
"final vertices after fitting: "<<endl;
2294 cout<<
" _mult_vtx " << i<<
" vertex "<<
_vertex_list[i][2]<<endl;
2299 for (
unsigned int k = 0;
k < nzvtx; ++
k){
2301 float zvtx = vtx_tracks[
k].z0;
2302 bool pass_dcaz=
false;
2305 if (new_vtx) izvtx=
iz;
2306 pass_dcaz = pass_dcaz || new_vtx;
2310 if (!pass_dcaz || izvtx >=
_vertex_list.size())
continue;
2325 for(
unsigned int i=0; i<vtx_tracks.size(); ++i) vtx_tracks[i].reset();
2327 _multi_vtx_tracks.clear();
2328 _multi_vtx_track_covars.clear();
2329 _multi_vtx_track_errors.clear();
2336 float phi,
float d,
float kappa,
float ,
float dzdl,
2337 Eigen::Matrix<float, 5, 5>
const& input,
2338 Eigen::Matrix<float, 6, 6>& output) {
2340 Eigen::Matrix<float, 6, 5> J = Eigen::Matrix<float, 6, 5>::Zero(6, 5);
2346 float nu = sqrt(kappa);
2347 float dk_dnu = 2 * nu;
2349 float cosphi = cos(phi);
2350 float sinphi = sin(phi);
2352 J(0, 0) = -sinphi *
d;
2354 J(1, 0) = d * cosphi;
2358 float pt = 0.003 * B / kappa;
2359 float dpt_dk = -0.003 * B / (kappa * kappa);
2361 J(3, 0) = -cosphi *
pt;
2362 J(3, 2) = -sinphi * dpt_dk * dk_dnu;
2363 J(4, 0) = -sinphi *
pt;
2364 J(4, 2) = cosphi * dpt_dk * dk_dnu;
2367 float alpha_half = pow(alpha, 0.5);
2368 float alpha_3_half = alpha * alpha_half;
2370 J(5, 2) = dpt_dk * dzdl * alpha_half * dk_dnu;
2371 J(5, 4) = pt * (alpha_half + dzdl * dzdl * alpha_3_half) * dk_dnu;
2373 output = J * input * (J.transpose());
2378 for (std::map<unsigned int,SimpleHit3D>::iterator
it=
hits_map.begin();
2390 for (
unsigned int itrack = 0; itrack <
_tracks.size(); ++itrack) {
2407 if (_phi < 0.) _phi += 2.*
M_PI;