63 #define LogDebug(exp) std::cout<<"DEBUG: " <<__FILE__<<": "<<__LINE__<<": "<< exp
64 #define LogError(exp) std::cout<<"ERROR: " <<__FILE__<<": "<<__LINE__<<": "<< exp
65 #define LogWarning(exp) std::cout<<"WARNING: "<<__FILE__<<": "<<__LINE__<<": "<< exp
73 using namespace Eigen;
77 unsigned int min_nlayers,
80 _t_output_io(nullptr),
83 _min_nlayers(min_nlayers),
85 _ca_min_nlayers(min_nlayers),
92 _use_max_kappa(
false),
93 fill_multi_zvtx(
true),
95 _max_kappa(numeric_limits<float>::
max()),
106 _ca_phi_cut(
M_PI/36.),
133 _bbc_vertexes(nullptr),
134 _clustermap(nullptr),
139 _hough_space(nullptr),
140 _hough_funcs(nullptr),
142 _ntp_zvtx_by_event(nullptr),
143 _ntp_zvtx_by_track(nullptr),
147 _kappa_d_phi(nullptr),
152 _layer_ilayer_map_all(),
167 separate_helicity(
false),
198 for (
unsigned int izoom =0; izoom<
nzooms ; ++izoom)
233 _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");
262 <<
"====================== PHPatternReco::InitRun() ======================"
264 cout <<
" Magnetic field set to: " <<
_mag_field <<
" Tesla" << endl;
265 cout <<
" Number of tracking layers: " <<
_nlayers << endl;
266 for (
unsigned int i = 0; i <
_nlayers; ++i) {
267 cout <<
" Tracking layer #" << i <<
" " <<
"radius = "
272 cout <<
" Minimum pT: " <<
_min_pt << endl;
273 cout <<
" Maximum DCA: " << boolalpha <<
_cut_on_dca << noboolalpha
276 cout <<
" Maximum DCA cut: " <<
_dcaxy_cut << endl;
278 cout <<
" Maximum DCAZ cut: " <<
_dcaz_cut << endl;
279 cout <<
" Momentum rescale factor: " <<
_pt_rescale << endl;
281 <<
"==========================================================================="
292 if (
Verbosity() > 0) cout <<
"PHPatternReco::process_event -- entered" << endl;
338 unsigned int nseq = 1;
339 unsigned int nattempt = 1;
345 if (iseq==nseq)
break;
358 for (
unsigned int iattempt =0; iattempt<nattempt; ++iattempt ){
360 cout<<iattempt <<
" th attempt "<<endl;
372 unsigned int zoomlevel = 0;
384 while(it_begin != it_end)
386 delete it_begin->second;
391 while(it_begin != it_end)
393 delete it_begin->second;
398 while(it_begin != it_end)
400 delete it_begin->second;
405 while(it_begin != it_end)
407 delete it_begin->second;
428 for (zoomlevel =1; zoomlevel<2; ++zoomlevel){
446 cout<<
"CellularAutomaton failed. "<<endl;
459 cout<<
"export output"<<endl;
500 void PHPatternReco::add_zoom(
unsigned int n_kappa,
unsigned int n_phi,
unsigned int n_d,
unsigned int n_dzdl,
unsigned int n_z0) {
501 std::vector<unsigned int> zoom {n_kappa, n_phi, n_d, n_dzdl, n_z0};
510 "PHCompositeNode",
"DST"));
512 cerr <<
PHWHERE <<
"DST Node missing, doing nothing." << endl;
525 cout <<
"SVTX node added" << endl;
530 "SvtxTrackMap",
"PHObject");
533 cout <<
"Svtx/SvtxTrackMap node added" << endl;
538 tb_node->
addNode(vertexes_node);
540 cout <<
"Svtx/SvtxVertexMap node added" << endl;
563 _radii.assign(_nlayers, 0.0);
564 map<float, int> radius_layer_map;
572 layerrange.first; layeriter != layerrange.second; ++layeriter) {
573 radius_layer_map.insert(
574 make_pair(layeriter->second->get_radius(),
575 layeriter->second->get_layer()));
581 laddergeos->get_begin_end();
583 layerrange.first; layeriter != layerrange.second; ++layeriter) {
584 radius_layer_map.insert(
585 make_pair(layeriter->second->get_radius(),
586 layeriter->second->get_layer()));
590 if (mvtxladdergeos) {
592 mvtxladdergeos->get_begin_end();
594 layerrange.first; layeriter != layerrange.second; ++layeriter) {
595 radius_layer_map.insert(
596 make_pair(layeriter->second->get_radius(),
597 layeriter->second->get_layer()));
606 for (map<float, int>::iterator iter = radius_layer_map.begin();
607 iter != radius_layer_map.end(); ++iter) {
624 for (; miter != begin_end.second; ++miter) {
646 laddergeos->get_begin_end();
648 for (; miter != begin_end.second; ++miter) {
666 if (mvtxladdergeos) {
668 mvtxladdergeos->get_begin_end();
670 for (; miter != begin_end.second; ++miter) {
692 map<int, float>::iterator mat_it;
714 _bbc_vertexes = findNode::getClass<BbcVertexMap>(topNode,
"BbcVertexMap");
724 _clustermap = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
727 cerr <<
PHWHERE <<
" ERROR: Can't find node TrkrClusterContainer" << endl;
730 _hitsets = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
733 std::cout <<
PHWHERE <<
"No hitset container on node tree. Bailing."
738 _trackmap = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMap");
740 cerr <<
PHWHERE <<
" ERROR: Can't find SvtxTrackMap." << endl;
745 _vertexmap = findNode::getClass<SvtxVertexMap>(topNode,
"SvtxVertexMap");
747 cerr <<
PHWHERE <<
" ERROR: Can't find SvtxVertexMap." << endl;
756 unsigned int clusid = 0;
757 unsigned int ilayer = 0;
760 for (
auto hitsetitr = hitsetrange.first;
761 hitsetitr != hitsetrange.second;
764 for(
auto clusIter = range.first; clusIter != range.second; ++clusIter ){
788 for (
int i = 0; i < 3; ++i) {
789 for (
int j = i; j < 3; ++j) {
796 cout <<
" inserting hit with layer " << ilayer <<
" into hits_map" << endl;
797 hits_map.insert(std::pair<unsigned int, SimpleHit3D>(hit3d.
get_id(),hit3d));
806 cout <<
"-------------------------------------------------------------------"
808 cout <<
"PHPatternReco::process_event has the following input clusters:"
811 cout <<
"n init clusters = " <<
hits_map.size() << endl;
813 for (std::map<unsigned int,SimpleHit3D>::iterator
it=
hits_map.begin();
822 cout <<
"-------------------------------------------------------------------"
867 std::map<int, SvtxVertex*> svtx_vertex_list;
876 svtx_vertex_list[vertex_id] = vertex;
881 for (
int i = 0; i < 3; ++i)
896 std::vector<SimpleHit3D> track_hits;
900 for (
unsigned int itrack = 0; itrack <
_tracks.size(); ++itrack) {
906 track_hits =
_tracks.at(itrack).hits;
908 for (
unsigned int ihit = 0;
ihit < track_hits.size(); ++
ihit) {
913 clusterkey = track_hits.at(
ihit).get_cluskey();
919 <<
": itrack: " << itrack
920 <<
": nhits: " << track_hits.size()
921 <<
": clusterkey: " << clusterkey
922 <<
": layer: " << layer
929 float kappa =
_tracks.at(itrack).kappa;
937 float x_center = cos(phi) * (d + 1 / kappa);
938 float y_center = sin(phi) * (d + 1 / kappa);
942 if ((track_hits[0].get_x()- x_center)*(track_hits[track_hits.size()-1].get_y()- y_center)
943 - (track_hits[0].get_y()- y_center)*(track_hits[track_hits.size()-1].get_x()- x_center)> 0)
952 pZ = pT * dzdl / sqrt(1.0 - dzdl * dzdl);
954 int ndf = 2 *
_tracks.at(itrack).hits.size() - 5;
957 track.
set_px(pT * cos(phi - helicity *
M_PI / 2));
958 track.
set_py(pT * sin(phi - helicity *
M_PI / 2));
977 for (
unsigned int row = 0; row < 6; ++row) {
978 for (
unsigned int col = 0;
col < 6; ++
col) {
983 unsigned int vid =9999;
984 float distance = 9999;
985 for (std::map<int,SvtxVertex*>::iterator
it=svtx_vertex_list.begin();
986 it!=svtx_vertex_list.end();
988 unsigned int iv =
it->second->get_id();
989 float zvtx =
it->second->get_position(2);
990 if (fabs(z0-zvtx) <distance){
992 distance = fabs(z0-zvtx);
997 if (fabs(z0 - distance)>0.05) vid = 9999;
998 if (
Verbosity() > 5) cout<<
"vertex_id "<<vid<<endl;
1005 track.
set_x(svtx_vertex_list[vid]->get_x() + d * cos(phi));
1006 track.
set_y(svtx_vertex_list[vid]->get_y() + d * sin(phi));
1007 track.
set_z(svtx_vertex_list[vid]->get_z() + z0);
1012 else svtx_vertex_list[vid]->insert_track(track.
get_id());
1015 cout <<
"track " << itrack <<
" quality = " << track.
get_quality() << endl;
1016 cout <<
"px = " << track.
get_px() <<
" py = " << track.
get_py() <<
" pz = " << track.
get_pz() << endl;
1068 for (
id=0;
id<
nd; ++
id) {
1081 bins_map.insert(std::pair<unsigned int, HelixHoughBin*>(
bin,_bin));
1089 cout<<
"bins_map.size " <<
bins_map.size()<<endl;
1096 std::vector<float> kappa_phi_d_ranges;
1103 kappa_phi_d_ranges.push_back(kappa_min);
1104 kappa_phi_d_ranges.push_back(kappa_max);
1105 kappa_phi_d_ranges.push_back(phi_min);
1106 kappa_phi_d_ranges.push_back(phi_max);
1107 kappa_phi_d_ranges.push_back(d_min);
1108 kappa_phi_d_ranges.push_back(d_max);
1115 for (std::map<unsigned int,SimpleHit3D>::iterator
it=
hits_map.begin();
1120 cluster_id =
it->first;
1124 auto hitused =
hits_used.find(cluster_id);
1126 used =
hits_used.find(cluster_id)->second;
1129 hitpos3d[0] = hit.
get_x();
1130 hitpos3d[1] = hit.
get_y();
1131 hitpos3d[2] = hit.
get_z();
1136 std::vector<float> z0_range;
1139 z0_range.push_back(z0_min);
1140 z0_range.push_back(z0_max);
1141 float dzdl_range[2] = {4,5};
1151 unsigned int dzdl_bin_range[2];
1154 if (dzdl_range[0] > dzdl_max)
continue;
1155 else if (dzdl_range[0] <dzdl_min) dzdl_bin_range[0]=0;
1156 if (dzdl_range[1] < dzdl_min)
continue;
1159 cout<<
"dzdl_min " <<dzdl_range[0]<<
" dzdl_max "<<dzdl_range[1]<<endl;
1161 for (
il =dzdl_bin_range[0];
il<dzdl_bin_range[1]+1; ++
il) {
1162 unsigned int zbins[5] = {0,0,0,
il,
iz};
1165 cout<<
"bin "<<
bin<<endl;
1170 bins_map.find(
bin)->second->add_cluster_ID(cluster_id);
1173 cout<<
"count "<<
bins_map.find(
bin)->second->get_count()<<endl;
1179 kappa_phi_d_ranges.clear();
1180 hitpos3d[0]=-999; hitpos3d[1]= -999; hitpos3d[2]=-999;
1183 cout<<
"total number of clusters "<<icluster<<endl;
1192 unsigned int max_i=-1;
1193 unsigned int max_count = 0;
1197 unsigned int bins[5]={0,0,0,j,i};
1199 unsigned int count =
bins_map.find(
bin)->second->get_count();
1200 if (count <1)
continue;
1201 if (count>max_count) {
1206 cout<<
"bin "<<
bin<<endl;
1207 cout<<
"z0 bin "<<i<<
" dzdl bin "<<j<<
" count "<< count<<endl;
1208 cout<<
"z0 bin "<<
bins_map.find(
bin)->second->get_z0_bin(0)<<
" dzdl bin "<<
bins_map.find(
bin)->second->get_dzdl_bin(0)<<endl;
1229 unsigned int bins[5]={0,0,0,
k,max_i};
1231 unsigned int count =
bins_map.find(
bin)->second->get_count();
1232 if (count <1)
continue;
1240 cout<<
"bins_map_sel.size at zoom "<<zoomlevel <<
" : (find_track_candidates_z_init) " <<
bins_map_sel.size()<<endl;
1272 hitpos3d[0] = hit.
get_x();
1273 hitpos3d[1] = hit.
get_y();
1274 hitpos3d[2] = hit.
get_z();
1279 std::vector<float> kappa_phi_d_ranges;
1286 kappa_phi_d_ranges.push_back(kappa_min);
1287 kappa_phi_d_ranges.push_back(kappa_max);
1288 kappa_phi_d_ranges.push_back(phi_min);
1289 kappa_phi_d_ranges.push_back(phi_max);
1290 kappa_phi_d_ranges.push_back(d_min);
1291 kappa_phi_d_ranges.push_back(d_max);
1302 std::vector<float> z0_range;
1305 z0_range.push_back(z0_min);
1306 z0_range.push_back(z0_max);
1307 float dzdl_range[2];
1311 unsigned int dzdl_bin_range[2];
1318 if (dzdl_range[0] > dzdl_max)
continue;
1319 else if (dzdl_range[0] <dzdl_min) dzdl_bin_range[0]=0;
1320 if (dzdl_range[1] < dzdl_min)
continue;
1325 for (
il =dzdl_bin_range[0];
il<dzdl_bin_range[1]+1; ++
il) {
1326 unsigned int zbins[5] = {0,0,0,
il,
iz};
1328 houghbin->
set_bin(zoomlevel,curbin);
1329 houghbin->
set_bins(zoomlevel,curbin);
1337 bins_map_cur.find(curgbin)->second->add_cluster_ID(cluster_id);
1344 cur_bin->
set_bin(zoomlevel,curbin);
1345 cur_bin->
set_bins(zoomlevel,curbin);
1347 bins_map_cur.insert(std::pair<unsigned int,HelixHoughBin*>(curgbin,cur_bin));
1355 unsigned int count =
bins_map_cur.find(curgbin)->second->get_count();
1362 hitpos3d[0]=-999; hitpos3d[1]= -999; hitpos3d[2]=-999;
1363 kappa_phi_d_ranges.clear();
1369 cout<<
"bins_map_cur.size at zoom "<<zoomlevel <<
" (vote_z) : " <<
bins_map_cur.size()<<endl;
1382 unsigned int count = houghbin->
get_count();
1383 if (count <1)
continue;
1400 cout<<
"il "<<
il<<
" iz "<<
iz<<
" count "<< count<<endl;
1408 cout<<
"bins_map_sel.size at zoom "<< zoomlevel<<
" (find_track_candidates_z) : " <<
bins_map_sel.size()<<endl;
1415 unsigned int ilfill= -1;
1416 unsigned int izfill= -1;
1422 unsigned int phi_bin_range[2];
1432 if (count==0){ ilfill =
il; izfill =
iz;}
1434 unsigned int ikprev =0;
1435 unsigned int ipprev =0;
1436 unsigned int idprev =0;
1440 idprev = houghbin->
get_d_bin(zoomlevel-1);
1443 cout<<
"bin "<<
bin<<endl;
1444 cout<<
"il " <<
il<<
" iz "<<
iz<<
" ikprev "<<ikprev<<
" ipprev "<<ipprev<<
" idprev "<<idprev<<
" count"<<houghbin->
get_count()<<endl;
1470 for (
id=0;
id<
nd; ++
id){
1471 std::vector<float> kappa_d_ranges;
1485 kappa_d_ranges.push_back(kappa_min);
1486 kappa_d_ranges.push_back(kappa_max);
1487 kappa_d_ranges.push_back(d_min);
1488 kappa_d_ranges.push_back(d_max);
1492 float phi_r_range[2];
1493 float phi_l_range[2];
1494 unsigned int phi_r_bin_range[2];
1495 unsigned int phi_l_bin_range[2];
1506 bool phi_r_out_of_range=
false;
1507 bool phi_l_out_of_range=
false;
1508 if (phi_r_range[0] > phi_max) phi_r_out_of_range=
true;
1509 else if (phi_r_range[0] <phi_min) phi_r_bin_range[0]=0;
1510 if (phi_r_range[1] < phi_min) phi_r_out_of_range=
true;
1513 if (phi_l_range[0] > phi_max) phi_r_out_of_range=
true;
1514 else if (phi_r_range[0] <phi_min) phi_r_bin_range[0]=0;
1515 if (phi_r_range[1] < phi_min) phi_r_out_of_range=
true;
1518 for (
int i = 0; i<2;++i){
1520 if ((i==0&&phi_r_out_of_range) || (i==1&&phi_l_out_of_range))
continue;
1523 phi_bin_range[0]=phi_r_bin_range[0];
1524 phi_bin_range[1]=phi_r_bin_range[1];
1527 phi_bin_range[0]=phi_l_bin_range[0];
1528 phi_bin_range[1]=phi_l_bin_range[1];
1533 for (
ip= phi_bin_range[0];
ip<phi_bin_range[1]+1;++
ip){
1536 houghbin->
set_bin(zoomlevel,curbin);
1537 houghbin->
set_bins(zoomlevel,curbin);
1550 cur_bin->
set_bin(zoomlevel,curbin);
1551 cur_bin->
set_bins(zoomlevel,curbin);
1553 bins_map_cur.insert(std::pair<unsigned int,HelixHoughBin*>(curgbin,cur_bin));
1564 kappa_d_phi_map.insert(std::pair<unsigned int,unsigned int>(curgbin,1));
1567 unsigned int kpbins[5]={
ik,
ip,0,
il,
iz};
1576 kappa_phi_map.insert(std::pair<unsigned int,unsigned int>(kpbin,1));
1579 unsigned int dpbins[5]={0,
ip,
id,
il,
iz};
1588 d_phi_map.insert(std::pair<unsigned int,unsigned int>(dpbin,1));
1592 if (
il ==ilfill &&
iz == izfill){
1599 cout <<
"il "<<
il<<
" iz "<<
iz<<endl;
1612 float phi_prev_range[2];
1613 float phi_next_range[2];
1614 unsigned int phi_bin_range[2];
1622 phi_prev_range[0]= phi_next_range[0];
1623 phi_prev_range[1]= phi_next_range[1];
1630 bool phi_out_of_range=
false;
1631 if (phi_range[0] > phi_max) phi_out_of_range=
true;
1632 else if (phi_range[0] <phi_min) phi_bin_range[0]=0;
1633 if (phi_range[1] < phi_min) phi_out_of_range=
true;
1636 if (phi_out_of_range)
continue;
1638 for (
ip= phi_bin_range[0];
ip<phi_bin_range[1]+1;++
ip){
1641 houghbin->
set_bin(zoomlevel,curbin);
1642 houghbin->
set_bins(zoomlevel,curbin);
1645 cout<<
" curgbin "<<curgbin<<endl;
1657 cur_bin->
set_bin(zoomlevel,curbin);
1658 cur_bin->
set_bins(zoomlevel,curbin);
1660 bins_map_cur.insert(std::pair<unsigned int,HelixHoughBin*>(curgbin,cur_bin));
1662 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;
1685 kappa_d_ranges.clear();
1691 cout<<
"bins_map_cur.size at zoom "<<zoomlevel <<
" (vote_xy) : " <<
bins_map_cur.size()<<endl;
1706 unsigned int count = houghbin->
get_count();
1707 if (count <1)
continue;
1716 cout<<
"bin "<<
bin<<endl;
1717 cout<<
"kappa bin "<<
ik<<
" phi bin "<<
ip<<
" d bin "<<
id <<
" il "<<
il<<
" iz "<<
iz<<
" count "<< count<<endl;
1725 cout<<
"inserting selected bin.. "<<endl;
1726 cout<<
"bin "<<
bin<<endl;
1727 cout<<
"kappa bin "<<
ik<<
" phi bin "<<
ip<<
" d bin "<<
id <<
" il "<<
il<<
" iz "<<
iz<<
" count "<< count<<endl;
1734 cout<<
"bins_map_prev.size at zoom "<<zoomlevel <<
" (find_track_candidates_xy) " <<
bins_map_prev.size()<<endl;
1744 unsigned int nclust = houghbin->
get_count();
1753 unsigned int var[3] = {1<<3,1<<4,(1<<3)+(1<<4)};
1754 unsigned int sign[3][4] ={{0<<3,1<<3,0,0},
1756 {(0<<3)+(0<<4),(0<<3)+(1<<4),(1<<3)+(0<<4),(1<<3)+(1<<4)}};
1757 for (
unsigned int lz =0; lz<3; ++lz){
1758 for (
unsigned int ss = 0; ss<4; ++ss ){
1759 if (lz<2 && ss>1)
continue;
1762 cout<<
"var "<<var[lz]<<
" ss "<<ss<<
" neighbor_global bin "<<neighborbin<<endl;
1764 if (globalbin==neighborbin)
continue;
1768 cout<<
"neighbor bin found"<<endl;
1770 unsigned int nclust_neighbor =
bins_map_sel.find(neighborbin)->second->get_count();
1771 if(nclust >= nclust_neighbor && nclust <= 1.1*
_nlayers) {
1773 cout<<
"merging "<<endl;
1776 iter !=
bins_map_sel.find(neighborbin)->second->end_clusters();
1790 cout<<
"bins_map_sel.size at zoom " <<zoomlevel<<
" (prune_z) : " <<
bins_map_sel.size()<<endl;
1796 cout<<
"bins_map_prev.size at zoom " <<zoomlevel<<
" (pre prune_xy) : " <<
bins_map_prev.size()<<endl;
1799 unsigned int nclust = houghbin->
get_count();
1802 cout<<
"global bin "<<globalbin<<endl;
1809 unsigned int var[6] = {1,1<<1,1<<2,(1)+(1<<1),(1)+(1<<2),(1<<1)+(1<<2)};
1810 unsigned int sign[6][4]={{0,1,0,0},
1813 {0,1<<1,1,(1)+(1<<1)},
1814 {0,1<<2,1, (1)+(1<<2)},
1815 {0,1<<2,1<<1,(1<<1)+(1<<2)}};
1816 for (
unsigned int kpd =0; kpd<6; ++kpd){
1817 for (
unsigned int ss = 0; ss<4; ++ss ){
1818 if (kpd<3 && ss>1)
continue;
1821 if (globalbin==neighborbin)
continue;
1825 unsigned int nclust_neighbor =
bins_map_prev.find(neighborbin)->second->get_count();
1826 if(nclust >= nclust_neighbor && nclust <= 1.1*
_nlayers) {
1828 cout<<
"merging "<<endl;
1831 iter !=
bins_map_prev.find(neighborbin)->second->end_clusters();
1845 cout<<
"bins_map_prev.size at zoom " <<zoomlevel<<
" (prune_xy) : " <<
bins_map_prev.size()<<endl;
1851 for (std::map<unsigned int,bool>::iterator
it =
hits_used.begin();
1859 unsigned int icluster = 0;
1871 track.
hits.push_back(hit);
1872 track.
hits[icluster].set_id(*iter);
1880 cout<<
"bins_map_prev selected, total bins : "<<
bins_map_prev.size()<<endl;
1896 auto search =
hits_map.find(*iter);
1900 track.
hits[icluster] = hit;
1901 track.
hits[icluster].set_id(*iter);
1911 new_tracks.push_back(track);
1915 cout<<
"Number of tracks added : (to be used for initial vertexing or track seeding)"<<new_tracks.size()<<endl;
1921 unsigned int layer0 = 0;
1922 unsigned int layer1 = 1;
1923 unsigned int layer2 = 2;
1930 std::vector<unsigned int> layer_sorted_0;
1931 std::vector<unsigned int> layer_sorted_1;
1932 std::vector<unsigned int> layer_sorted_2;
1934 for (std::map<unsigned int,SimpleHit3D>::iterator
it=
hits_map.begin();
1948 unsigned int hitid = hit.
get_id();
1949 if (hitid != hitid)
continue;
1951 if (layer == layer0 ) layer_sorted_0.push_back(hitid);
1952 else if (layer == layer1) layer_sorted_1.push_back(hitid);
1953 else if (layer == layer2) layer_sorted_2.push_back(hitid);
1957 cout<<
"layer 0 " <<layer_sorted_0.size()<<endl;
1958 cout<<
"layer 1 " <<layer_sorted_1.size()<<endl;
1959 cout<<
"layer 2 " <<layer_sorted_2.size()<<endl;
1962 std::set<unsigned int> clusters;
1963 bool fill_track =
false;
1964 unsigned int cluster_layer0 = 0;
1965 for (std::vector<unsigned int>::iterator it0 = layer_sorted_0.begin() ; it0 != layer_sorted_0.end(); ++it0)
1974 auto search0 =
hits_map.find(*it0);
1975 if (search0 ==
hits_map.end())
continue;
1977 float phi_h0 = atan2(hit0.
get_y(),hit0.
get_x());
1979 float z_h0 = hit0.
get_z();
1983 clusters.insert(*it0);
1985 for (std::vector<unsigned int>::iterator it1 = layer_sorted_1.begin(); it1 != layer_sorted_1.end(); ++it1 )
1988 float phi_h1 = atan2(hit1.
get_y(),hit1.
get_x());
1990 float z_h1 = hit1.
get_z();
1992 float phi_h0h1 = phi_h1-phi_h0;
1993 if (phi_h1< M_PI/2. && phi_h0 > 3*
M_PI/2.) phi_h0h1 += 2.*
M_PI;
1994 else if (phi_h1> 3*
M_PI/2. && phi_h0 <
M_PI/2.) phi_h0h1 -= 2.*
M_PI;
1999 for(std::vector<unsigned int>::iterator it2 = layer_sorted_2.begin(); it2 != layer_sorted_2.end(); ++it2)
2002 float phi_h2 = atan2(hit2.
get_y(),hit2.
get_x());
2004 float phi_h1h2 = phi_h2-phi_h1;
2005 if (phi_h2< M_PI/2. && phi_h1 > 3*
M_PI/2.) phi_h1h2 += 2.*
M_PI;
2006 else if (phi_h2> 3*
M_PI/2. && phi_h1 <
M_PI/2.) phi_h1h2 -= 2.*
M_PI;
2009 float z_h2 = hit2.
get_z();
2010 if (fabs(phi_h1h2) <
_ca_phi_cut && (z_h1-z_h0)*(z_h2-z_h1)>0 && fabs(z_h2-z_h1)<
_z_cut ){
2012 clusters.insert(*it1);
2013 clusters.insert(*it2);
2021 unsigned int nclusters =0 ;
2023 for (std::set<unsigned int>::iterator
it=clusters.begin();
it!=clusters.end(); ++
it ){
2026 track.
hits[nclusters] = hit;
2027 track.
hits[nclusters].set_id(*
it);
2030 new_tracks.push_back(track);
2035 cout<<
"number of triplets "<<new_tracks.size()<<endl;
2052 auto hitused =
hits_used.find(cluster_id);
2054 hits_used.find(cluster_id)->second =
true;
2063 cout<<
"Entering cellular autumaton : processing "<< candidate_tracks.size()<<
" tracks. "<<endl;
2087 for(
unsigned int i=0; i<candidate_tracks.size(); ++i) candidate_tracks[i].reset();
2088 candidate_tracks.clear();
2149 cout<<
"tracks size " <<
_tracks.size()<<endl;
2151 if (
_tracks.empty())
return -1;
2152 std::vector<SimpleTrack3D> vtx_tracks;
2156 unsigned int nzvtx = vtx_tracks.size();
2159 unsigned int nzbins= 2*1280;
2161 float zvalues[2*1280];
2162 unsigned int zcounts[2*1280];
2163 for (
unsigned int i = 0; i<nzbins; ++i)
2165 zvalues[i] =
_min_z0 + (i+0.5)*binsize;
2169 for (
unsigned int i = 0; i < nzvtx; ++i)
2171 double zvtx = vtx_tracks[i].z0;
2172 if (zvtx != zvtx || zvtx<_min_z0 || zvtx>
_max_z0)
continue;
2173 unsigned int zbin = (zvtx-
_min_z0)/binsize;
2174 cout<<
"zvtx " <<zvtx <<
" zbin "<< zbin <<endl;
2178 std::vector<float> zvertices;
2179 std::vector<float> zmeans;
2180 std::vector<float> zsums;
2181 std::vector<float> zsigmas;
2182 std::vector<unsigned int> nzvtxes;
2184 for (
unsigned int j=2; j<(nzbins-3); ++j)
2186 cout<<
"z countz "<<j<<
" "<<zcounts[j]<<endl;
2189 bool twobinspeak =
_mult_twobins*(zcounts[j-1]+zcounts[j-2])<(zcounts[j]+zcounts[j+1])
2190 &&
_mult_twobins*(zcounts[j+2]+zcounts[j+3])< (zcounts[j]+zcounts[j+1]);
2196 float zvertex=-999.;
2200 if (onebinpeak) zvertex = zvalues[j];
2201 if (twobinspeak) zvertex = (zvalues[j]*zcounts[j] + zvalues[j+1]*zcounts[j+1])/(zcounts[j]+zcounts[j+1]);
2202 zvertices.push_back(zvertex);
2203 nzvtxes.push_back(0);
2204 zmeans.push_back(0.0);
2205 zsums.push_back(0.0);
2206 zsigmas.push_back(0.0);
2207 cout<<
"vertex "<< zvertex<<endl;
2211 cout<<
"number of candidate z-vertices : "<< zvertices.size()<<endl;
2212 for (
unsigned int j = 0; j <zvertices.size(); ++j)
2213 cout<<
"primary vertex position : "<<zvertices[j]<<endl;
2215 unsigned int izvtx= 999;
2216 for (
unsigned int i = 0; i < nzvtx; ++i) {
2217 double zvtx = vtx_tracks[i].z0;
2218 if (zvtx != zvtx)
continue;
2219 bool pass_dcaz=
false;
2220 for (
unsigned int k = 0;
k<zvertices.size(); ++
k) {
2221 bool new_vtx = fabs(zvtx-zvertices[
k]) <
_dcaz_cut;
2222 if (new_vtx) izvtx=
k;
2223 pass_dcaz = pass_dcaz || new_vtx;
2226 if (!pass_dcaz || izvtx>99)
continue;
2228 zsums[izvtx] += zvtx;
2231 for (
unsigned int iz = 0;
iz<zvertices.size();++
iz ){
2232 if (nzvtxes[
iz]==0)
continue;
2233 zmeans[
iz] = zsums[
iz]/nzvtxes[
iz];
2234 cout<<
"zmean for vertex "<<
iz <<
" = "<<zmeans[
iz]<<endl;
2235 zvertices[
iz] = zmeans[
iz];
2240 ntp_data[1] = zmeans[
iz];
2245 for (
unsigned int j = 0; j < nzvtx; ++j){
2246 double zvtx = vtx_tracks[j].z0;
2247 if (zvtx != zvtx)
continue;
2253 ntp_data[1] = vtx_tracks[j].z0;
2259 bool pass_dcaz=
false;
2260 for (
unsigned int k = 0;
k<zvertices.size(); ++
k) {
2261 bool new_vtx = fabs(zvtx-zvertices[
k]) <
_dcaz_cut;
2262 if (new_vtx) izvtx=
k;
2263 pass_dcaz = pass_dcaz || new_vtx;
2266 if (!pass_dcaz || izvtx>99)
continue;
2267 zsums[izvtx] += pow(fabs(zmeans[izvtx] - vtx_tracks[j].
z0),2);
2270 std::vector<float> _multi_vtx;
2271 std::vector<std::vector<SimpleTrack3D> > _multi_vtx_tracks;
2272 std::vector<std::vector<double> > _multi_vtx_track_errors;
2273 std::vector<std::vector<Eigen::Matrix<float, 5, 5> > > _multi_vtx_track_covars;
2276 for (
unsigned int i = 0; i < zvertices.size(); ++i)
2278 if (nzvtxes[i]==0)
continue;
2279 _multi_vtx.push_back(zvertices[i]);
2280 std::vector<SimpleTrack3D> one_vtx_tracks;
2281 _multi_vtx_tracks.push_back(one_vtx_tracks);
2282 std::vector<double> one_track_errors;
2283 _multi_vtx_track_errors.push_back(one_track_errors);
2284 std::vector<Eigen::Matrix<float, 5, 5> > one_track_covars;
2285 _multi_vtx_track_covars.push_back(one_track_covars);
2286 zsigmas[i] = sqrt(zsums[i]/(nzvtxes[i]-1));
2287 cout<<
"zsigma for vertex "<<i <<
" = "<<zsigmas[i]<<endl;
2290 unsigned int nzvtx_final = 0;
2291 for (
unsigned int k = 0;
k < nzvtx; ++
k){
2293 float zvtx = vtx_tracks[
k].z0;
2294 bool pass_dcaz=
false;
2295 for (
unsigned int iz = 0;
iz<_multi_vtx.size(); ++
iz) {
2296 bool new_vtx = fabs(zvtx-_multi_vtx[
iz]) <
_dcaz_cut;
2297 if (new_vtx) izvtx=
iz;
2298 pass_dcaz = pass_dcaz || new_vtx;
2301 if (!pass_dcaz || izvtx >= _multi_vtx.size())
continue;
2306 _multi_vtx_tracks[izvtx].push_back(vtx_tracks[
k]);
2308 _multi_vtx_track_errors[izvtx].push_back(
_track_states[k].chi2);
2312 cout<<
"start fitting vertex.. "<<endl;
2313 for (
unsigned int i = 0; i<_multi_vtx.size(); ++i)
2315 if (_multi_vtx_tracks[i].size()==0)
continue;
2318 cout <<
" seed track vertex pre-fit: "
2329 cout<<
"number of fitted tracks for vertex "<<i<<
" : "<<
n_vtx_tracks <<endl;
2332 cout <<
" seed track vertex post-fit: "
2343 cout<<
"final vertices : "<<endl;
2349 for (
unsigned int k = 0;
k < nzvtx; ++
k){
2351 float zvtx = vtx_tracks[
k].z0;
2352 bool pass_dcaz=
false;
2355 if (new_vtx) izvtx=
iz;
2356 pass_dcaz = pass_dcaz || new_vtx;
2360 if (!pass_dcaz || izvtx >=
_vertex_list.size())
continue;
2374 for(
unsigned int i=0; i<vtx_tracks.size(); ++i) vtx_tracks[i].reset();
2376 _multi_vtx_tracks.clear();
2377 _multi_vtx_track_covars.clear();
2378 _multi_vtx_track_errors.clear();
2385 float phi,
float d,
float kappa,
float z0,
float dzdl,
2386 Eigen::Matrix<float, 5, 5>
const& input,
2387 Eigen::Matrix<float, 6, 6>& output) {
2389 Eigen::Matrix<float, 6, 5> J = Eigen::Matrix<float, 6, 5>::Zero(6, 5);
2395 float nu = sqrt(kappa);
2396 float dk_dnu = 2 * nu;
2398 float cosphi = cos(phi);
2399 float sinphi = sin(phi);
2401 J(0, 0) = -sinphi *
d;
2403 J(1, 0) = d * cosphi;
2407 float pt = 0.003 * B / kappa;
2408 float dpt_dk = -0.003 * B / (kappa * kappa);
2410 J(3, 0) = -cosphi *
pt;
2411 J(3, 2) = -sinphi * dpt_dk * dk_dnu;
2412 J(4, 0) = -sinphi *
pt;
2413 J(4, 2) = cosphi * dpt_dk * dk_dnu;
2416 float alpha_half = pow(alpha, 0.5);
2417 float alpha_3_half = alpha * alpha_half;
2419 J(5, 2) = dpt_dk * dzdl * alpha_half * dk_dnu;
2420 J(5, 4) = pt * (alpha_half + dzdl * dzdl * alpha_3_half) * dk_dnu;
2422 output = J * input * (J.transpose());
2427 for (std::map<unsigned int,SimpleHit3D>::iterator
it=
hits_map.begin();
2439 for (
unsigned int itrack = 0; itrack <
_tracks.size(); ++itrack) {
2456 if (_phi < 0.) _phi += 2.*
M_PI;