43 #include <boost/geometry.hpp>
44 #include <boost/geometry/geometries/box.hpp>
45 #include <boost/geometry/geometries/point.hpp>
46 #include <boost/geometry/index/rtree.hpp>
60 typedef bg::model::point<float, 3, bg::cs::cartesian>
point;
61 typedef bg::model::box<point>
box;
62 typedef std::pair<point, TrkrDefs::cluskey>
pointKey;
64 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
65 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
66 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
70 namespace bg = boost::geometry;
71 namespace bgi = boost::geometry::index;
77 unsigned int nlayers_maps,
78 unsigned int nlayers_intt,
79 unsigned int nlayers_tpc,
80 unsigned int start_layer)
82 , _g4clusters(nullptr)
84 , _g4vertexes(nullptr)
85 , _cluster_map(nullptr)
86 , _svtxhitsmap(nullptr)
87 , _hit_used_map(nullptr)
88 , _hit_used_map_size(0)
95 , _nlayers_maps(nlayers_maps)
96 , _nlayers_intt(nlayers_intt)
97 , _nlayers_tpc(nlayers_tpc)
98 , _start_layer(start_layer)
115 map<float, int> radius_layer_map;
126 layeriter != layerrange.second; ++layeriter)
128 radius_layer_map.insert(
129 make_pair(layeriter->second->get_radius(),
130 layeriter->second->get_layer()));
137 laddergeos->get_begin_end();
140 layeriter != layerrange.second; ++layeriter)
142 radius_layer_map.insert(
143 make_pair(layeriter->second->get_radius(),
144 layeriter->second->get_layer()));
151 mapsladdergeos->get_begin_end();
154 layeriter != layerrange.second; ++layeriter)
156 radius_layer_map.insert(
157 make_pair(layeriter->second->get_radius(),
158 layeriter->second->get_layer()));
161 for (map<float, int>::iterator iter = radius_layer_map.begin();
162 iter != radius_layer_map.end(); ++iter)
178 for (; miter != begin_end.second; miter++)
195 laddergeos->get_begin_end();
197 for (; miter != begin_end.second; miter++)
213 mapsladdergeos->get_begin_end();
215 for (; miter != begin_end.second; miter++)
249 const vector<TrkrCluster *> &
v2)
251 std::vector<TrkrCluster *>
v3;
252 vector<TrkrCluster *> va(v1);
253 vector<TrkrCluster *> vb(v2);
256 std::set_intersection(va.begin(), va.end(),
257 vb.begin(), vb.end(),
260 if (v3.size() < 3)
return false;
279 vector<TrkrCluster *>
yunion(
const vector<TrkrCluster *> &
v1,
const vector<TrkrCluster *> &
v2)
281 std::vector<TrkrCluster *>
v3;
286 std::set_union(v1.begin(), v1.end(),
287 v2.begin(), v2.end(),
293 const vector<TrkrCluster *> &
smaller)
295 for (
unsigned int i = 0; i < smaller.size(); i++)
297 if (!binary_search(larger.begin(), larger.end(), smaller.at(i),
comparing))
return false;
303 const vector<TrkrCluster *> &
v2)
305 vector<TrkrCluster *> temp(
yunion(v1, v2));
322 double predictedx = xx[0] - sin(xx[2]) / xx[4] + cos(xx[4] / tan(xx[3]) * qz + xx[2] -
M_PI / 2) / xx[4];
323 double predictedy = xx[1] + cos(xx[2]) / xx[4] + sin(xx[4] / tan(xx[3]) * qz + xx[2] -
M_PI / 2) / xx[4];
324 double xresid = tc->
getX() - predictedx;
325 double yresid = tc->
getY() - predictedy;
326 double zresid = tc->
getZ() - xx[i + 5];
328 double contribution = xresid * xresid / tc->
getError(0, 0) + yresid * yresid / tc->
getError(1, 1) + zresid * zresid / tc->
getError(2, 2);
330 chi2 += contribution;
333 double predictedx = xx[0] - sin(xx[2]) / xx[4] + cos(xx[4] / tan(xx[3]) * qz + xx[2] -
M_PI / 2) / xx[4];
334 double predictedy = xx[1] + cos(xx[2]) / xx[4] + sin(xx[4] / tan(xx[3]) * qz + xx[2] -
M_PI / 2) / xx[4];
341 chi2 += contribution;
351 _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
354 cerr <<
PHWHERE <<
" ERROR: Can't find node TRKR_CLUSTER" << endl;
362 double s = phi1 + phi2;
373 double d = phi1 - phi2;
382 void PHRTreeSeeding::QueryTree(
const bgi::rtree<
pointKey, bgi::quadratic<16>> &rtree,
double phimin,
double etamin,
double lmin,
double phimax,
double etamax,
double lmax, std::vector<pointKey> &returned_values)
384 rtree.query(bgi::intersects(
box(
point(phimin, etamin, lmin),
point(phimax, etamax, lmax))), std::back_inserter(returned_values));
385 if (phimin < 0) rtree.query(bgi::intersects(
box(
point(2 *
M_PI + phimin, etamin, lmin),
point(2 *
M_PI, etamax, lmax))), std::back_inserter(returned_values));
386 if (phimax > 2 *
M_PI) rtree.query(bgi::intersects(
box(
point(0, etamin, lmin),
point(phimax - 2 *
M_PI, etamax, lmax))), std::back_inserter(returned_values));
416 return xx[0] + xx[1] + xx[2] + xx[3] + xx[4];
430 double predictedx = xx[0] - sin(xx[2]) / xx[4] + cos(xx[4] / tan(xx[3]) * qz + xx[2] -
M_PI / 2) / xx[4];
431 double predictedy = xx[1] + cos(xx[2]) / xx[4] + sin(xx[4] / tan(xx[3]) * qz + xx[2] -
M_PI / 2) / xx[4];
432 double qx = tc->
getX() - predictedx;
433 double qy = tc->
getY() - predictedy;
435 chi2 += (qd * qx * qx - (qb + qc) * qx * qy + qa * qy * qy) / (qa * qd - qb * qc);
448 for (
int j = 0; j < 60; j++) nlayer[j] = 0;
451 for (
auto hitsetitr = hitsetrange.first;
452 hitsetitr != hitsetrange.second;
455 for(
auto clusIter = range.first; clusIter != range.second; ++clusIter ){
459 if (layer < 39)
continue;
463 double clus_phi = vec.Phi();
464 clus_phi -= 2 *
M_PI * floor(clus_phi / (2 *
M_PI));
465 double clus_eta = vec.Eta();
466 double clus_l =
layer;
468 vector<pointKey> testduplicate;
469 QueryTree(
_rtree, clus_phi - 0.00001, clus_eta - 0.00001, layer - 0.5, clus_phi + 0.00001, clus_eta + 0.00001, layer + 0.5, testduplicate);
470 if (!testduplicate.empty())
477 _rtree.insert(std::make_pair(
point(clus_phi, clus_eta, clus_l), ckey));
482 std::cout <<
"number of duplicates : " << n_dupli << std::endl;
487 TFile fpara(
"rtree_para.root",
"RECREATE");
488 TNtuple *NT =
new TNtuple(
"NT",
"NT",
"pt:dpt:z:dz:phi:dphi:c:dc:nhit");
511 int numberofseeds = 0;
512 cout <<
" entries in tree: " <<
_rtree.size() << endl;
514 for (
unsigned int iteration = 0; iteration < 1; ++iteration)
517 vector<pointKey> StartLayerClusters;
518 _rtree.query(bgi::intersects(
521 std::back_inserter(StartLayerClusters));
523 for (vector<pointKey>::iterator StartCluster = StartLayerClusters.begin(); StartCluster != StartLayerClusters.end(); ++StartCluster)
525 double StartPhi = StartCluster->first.get<0>();
526 double StartEta = StartCluster->first.get<1>();
528 vector<pointKey> SecondLayerClusters;
536 SecondLayerClusters);
537 cout <<
" entries in second layer: " << SecondLayerClusters.size() << endl;
539 for (vector<pointKey>::iterator SecondCluster = SecondLayerClusters.begin(); SecondCluster != SecondLayerClusters.end(); ++SecondCluster)
541 double currentphi = SecondCluster->first.get<0>();
542 double currenteta = SecondCluster->first.get<1>();
552 vector<double> curvatureestimates;
553 curvatureestimates.push_back(copysign(2 / sqrt(ther * ther + 1 / dphidr / dphidr), dphidr));
554 vector<double> phi_zigzag;
555 phi_zigzag.push_back(dphidr);
556 vector<double> z_zigzag;
557 z_zigzag.push_back(detadr);
559 vector<TrkrDefs::cluskey> cluskeys;
560 cluskeys.push_back(StartCluster->second);
561 cluskeys.push_back(SecondCluster->second);
562 cout <<
" phi 1: " << StartPhi
563 <<
" phi 2: " << currentphi
564 <<
" dphi: " << dphidr
565 <<
" eta 1: " << StartEta
566 <<
" eta 2: " << currenteta
567 <<
" deta: " << detadr
572 vector<pointKey> newlayer_clusters;
576 <<
" etamin " << currenteta -
etast
577 <<
" etamax " << currenteta +
etast
588 if (newlayer_clusters.empty())
591 if (failures > 2)
break;
595 double xinrecord = 100.0;
596 pointKey *xinkey = &*newlayer_clusters.begin();
597 for (std::vector<pointKey>::iterator
it = newlayer_clusters.begin();
it != newlayer_clusters.end(); ++
it)
601 cout <<
" nuphi: " <<
it->first.get<0>()
602 <<
" nueta: " <<
it->first.get<1>()
604 <<
" lay: " << newlayer
605 <<
" dl: " << lastgoodlayer - newlayer
609 if (dist < xinrecord)
617 detadr = (currenteta - xinkey->first.get<1>()) / (
_radii_all[lastgoodlayer] -
_radii_all[newlayer]);
620 curvatureestimates.push_back(copysign(2 / sqrt(ther * ther + 1 / dphidr / dphidr), dphidr));
622 phi_zigzag.push_back(dphidr);
623 z_zigzag.push_back(detadr);
625 cluskeys.push_back(xinkey->second);
627 currentphi = xinkey->first.get<0>();
628 currenteta = (currenteta + xinkey->first.get<1>()) / 2;
630 lastgoodlayer = newlayer;
633 if (failures > 2)
continue;
635 double phi_sum = std::accumulate(phi_zigzag.begin(), phi_zigzag.end(), 0.0);
636 double phi_mean = phi_sum / phi_zigzag.size();
638 std::vector<double> phi_diff(phi_zigzag.size());
639 std::transform(phi_zigzag.begin(), phi_zigzag.end(), phi_diff.begin(),
640 std::bind2nd(std::minus<double>(), phi_mean));
641 double phi_sq_sum = std::inner_product(phi_diff.begin(), phi_diff.end(), phi_diff.begin(), 0.0);
642 double phi_stdev = std::sqrt(phi_sq_sum / (phi_zigzag.size() - 1));
644 double z_sum = std::accumulate(z_zigzag.begin(), z_zigzag.end(), 0.0);
645 double z_mean = z_sum / z_zigzag.size();
647 std::vector<double> z_diff(z_zigzag.size());
649 std::bind2nd(std::minus<double>(), z_mean));
650 double z_sq_sum = std::inner_product(z_diff.begin(), z_diff.end(), z_diff.begin(), 0.0);
651 double z_stdev = std::sqrt(z_sq_sum / (z_zigzag.size() - 1));
653 double curv_sum = std::accumulate(curvatureestimates.begin(), curvatureestimates.end(), 0.0);
654 double curv_mean = curv_sum / curvatureestimates.size();
656 std::vector<double> curv_diff(curvatureestimates.size());
657 std::transform(curvatureestimates.begin(), curvatureestimates.end(), curv_diff.begin(),
658 std::bind2nd(std::minus<double>(), curv_mean));
659 double curv_sq_sum = std::inner_product(curv_diff.begin(), curv_diff.end(), curv_diff.begin(), 0.0);
660 double curv_stdev = std::sqrt(curv_sq_sum / (curvatureestimates.size() - 1));
662 const double BQ = 0.01 * 1.4 * 0.299792458;
663 double pt = BQ /
abs(curv_mean);
664 double pterror = BQ * curv_stdev / (curv_mean * curv_mean);
667 NT->Fill(pt, pterror, z_mean, z_stdev, phi_mean, phi_stdev, curv_mean, curv_stdev, cluskeys.size());
669 track.
set_id(numberofseeds);
671 for (
unsigned int j = 0; j < cluskeys.size(); j++)
675 track.
set_ndf(2 * cluskeys.size() - 5);
676 short int helicity = 1;
677 if (StartPhi * curv_mean < 0) helicity -= 2;
684 track.
set_px(pt * cos(StartPhi));
685 track.
set_py(pt * sin(StartPhi));
686 track.
set_pz(pt / tan(2 * atan(exp(-StartEta))));
693 track.
set_error(3, 3, pterror * pterror * cos(StartPhi) * cos(StartPhi));
694 track.
set_error(4, 4, pterror * pterror * sin(StartPhi) * sin(StartPhi));
695 track.
set_error(5, 5, pterror * pterror / tan(2 * atan(exp(-StartEta))) / tan(2 * atan(exp(-StartEta))));
703 cout <<
"number of seeds " << numberofseeds << endl;
719 cout <<
"Called Setup" << endl;
720 cout <<
"topNode:" << topNode << endl;
731 cout <<
"Called End " << endl;