Warning, /coresoftware/offline/packages/trackreco/PHHybridSeeding.cc.outdated is written in an unsupported language. File is not indexed.
0001 /*!
0002 * \file PHHybridSeeding.cc
0003 * \brief Track Seeding using STAR "CA" algorithm and ALICE simplified Kalman filter
0004 * \detail
0005 * \author Michael Peters & Christof Roland
0006 */
0007
0008 #include "PHHybridSeeding.h"
0009 #include "GPUTPCTrackLinearisation.h"
0010 #include "GPUTPCTrackParam.h"
0011
0012 #include "../PHTpcTracker/externals/kdfinder.hpp"
0013 #include "../PHTpcTracker/PHTpcTrackerUtil.h"
0014
0015 // trackbase_historic includes
0016 #include <trackbase_historic/SvtxTrackMap.h>
0017 #include <trackbase_historic/SvtxTrack_v3.h>
0018 #include <trackbase_historic/SvtxTrack.h>
0019 #include <globalvertex/SvtxVertex.h>
0020 #include <globalvertex/SvtxVertexMap.h>
0021
0022 #include <trackbase/TrkrCluster.h> // for TrkrCluster
0023 #include <trackbase/TrkrClusterContainer.h>
0024 #include <trackbase/TrkrDefs.h> // for getLayer, clu...
0025
0026 // sPHENIX Geant4 includes
0027 #include <g4detectors/PHG4TpcCylinderGeom.h>
0028 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0029 #include <g4detectors/PHG4CylinderGeom.h>
0030 #include <g4detectors/PHG4CylinderGeomContainer.h>
0031
0032 // sPHENIX includes
0033 #include <fun4all/Fun4AllReturnCodes.h>
0034
0035 #include <phool/PHTimer.h> // for PHTimer
0036 #include <phool/getClass.h>
0037 #include <phool/phool.h> // for PHWHERE
0038
0039 //ROOT includes
0040 #include <TVector3.h> // for TVector3
0041 #include <TFile.h>
0042 #include <TNtuple.h>
0043
0044 #include <Eigen/Core>
0045 #include <Eigen/Dense>
0046
0047 #include <algorithm>
0048 #include <cmath>
0049 #include <iostream>
0050 #include <numeric>
0051 #include <utility> // for pair, make_pair
0052 #include <vector>
0053 #include <algorithm> // for find
0054 #include <unordered_set>
0055
0056 // forward declarations
0057 class PHCompositeNode;
0058
0059 #define _DEBUG_
0060
0061 #if defined(_DEBUG_)
0062 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
0063 #else
0064 #define LogDebug(exp) (void)0
0065 #endif
0066
0067 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
0068 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
0069
0070 //end
0071
0072 typedef std::vector<TrkrDefs::cluskey> keylist;
0073 using namespace std;
0074
0075 PHHybridSeeding::PHHybridSeeding(
0076 const string &name,
0077 double maxSinPhi,
0078 double fieldDir,
0079 double search_radius1,
0080 double search_angle1,
0081 size_t min_track_size1,
0082 double search_radius2,
0083 double search_angle2,
0084 size_t min_track_size2,
0085 size_t nthreads
0086 )
0087 : PHTrackSeeding(name)
0088 , _max_sin_phi(maxSinPhi)
0089 , _fieldDir(fieldDir)
0090 , _search_radius1(search_radius1)
0091 , _search_angle1(search_angle1)
0092 , _min_track_size1(min_track_size1)
0093 , _search_radius2(search_radius2)
0094 , _search_angle2(search_angle2)
0095 , _min_track_size2(min_track_size2)
0096 , _nthreads(nthreads)
0097 {
0098 }
0099
0100 int PHHybridSeeding::InitializeGeometry(PHCompositeNode *topNode)
0101 {
0102 PHG4TpcCylinderGeomContainer *cellgeos = findNode::getClass<
0103 PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0104 PHG4CylinderGeomContainer *laddergeos = findNode::getClass<
0105 PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0106 PHG4CylinderGeomContainer *mapsladdergeos = findNode::getClass<
0107 PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
0108
0109 //_nlayers_seeding = _seeding_layer.size();
0110 //_radii.assign(_nlayers_seeding, 0.0);
0111 map<double, int> radius_layer_map;
0112
0113 _radii_all.assign(60, 0.0);
0114 _layer_ilayer_map.clear();
0115 _layer_ilayer_map_all.clear();
0116 if (cellgeos)
0117 {
0118 PHG4TpcCylinderGeomContainer::ConstRange layerrange =
0119 cellgeos->get_begin_end();
0120 for (PHG4TpcCylinderGeomContainer::ConstIterator layeriter =
0121 layerrange.first;
0122 layeriter != layerrange.second; ++layeriter)
0123 {
0124 radius_layer_map.insert(
0125 make_pair(layeriter->second->get_radius(),
0126 layeriter->second->get_layer()));
0127 }
0128 }
0129
0130 if (laddergeos)
0131 {
0132 PHG4CylinderGeomContainer::ConstRange layerrange =
0133 laddergeos->get_begin_end();
0134 for (PHG4CylinderGeomContainer::ConstIterator layeriter =
0135 layerrange.first;
0136 layeriter != layerrange.second; ++layeriter)
0137 {
0138 radius_layer_map.insert(
0139 make_pair(layeriter->second->get_radius(),
0140 layeriter->second->get_layer()));
0141 }
0142 }
0143
0144 if (mapsladdergeos)
0145 {
0146 PHG4CylinderGeomContainer::ConstRange layerrange =
0147 mapsladdergeos->get_begin_end();
0148 for (PHG4CylinderGeomContainer::ConstIterator layeriter =
0149 layerrange.first;
0150 layeriter != layerrange.second; ++layeriter)
0151 {
0152 radius_layer_map.insert(
0153 make_pair(layeriter->second->get_radius(),
0154 layeriter->second->get_layer()));
0155 }
0156 }
0157 for (map<double, int>::iterator iter = radius_layer_map.begin();
0158 iter != radius_layer_map.end(); ++iter)
0159 {
0160 _layer_ilayer_map_all.insert(make_pair(iter->second, _layer_ilayer_map_all.size()));
0161
0162 /*if (std::find(_seeding_layer.begin(), _seeding_layer.end(),
0163 iter->second) != _seeding_layer.end())
0164 {
0165 _layer_ilayer_map.insert(make_pair(iter->second, ilayer));
0166 ++ilayer;
0167 }*/
0168 }
0169 if (cellgeos)
0170 {
0171 PHG4TpcCylinderGeomContainer::ConstRange begin_end =
0172 cellgeos->get_begin_end();
0173 PHG4TpcCylinderGeomContainer::ConstIterator miter = begin_end.first;
0174 for (; miter != begin_end.second; ++miter)
0175 {
0176 PHG4TpcCylinderGeom *geo = miter->second;
0177 _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
0178 geo->get_radius() + 0.5 * geo->get_thickness();
0179
0180 /*if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
0181 {
0182 _radii[_layer_ilayer_map[geo->get_layer()]] =
0183 geo->get_radius();
0184 }*/
0185 }
0186 }
0187
0188 if (laddergeos)
0189 {
0190 PHG4CylinderGeomContainer::ConstRange begin_end =
0191 laddergeos->get_begin_end();
0192 PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
0193 for (; miter != begin_end.second; ++miter)
0194 {
0195 PHG4CylinderGeom *geo = miter->second;
0196 _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
0197 geo->get_radius() + 0.5 * geo->get_thickness();
0198
0199 /*if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
0200 {
0201 _radii[_layer_ilayer_map[geo->get_layer()]] = geo->get_radius();
0202 }*/
0203 }
0204 }
0205
0206 if (mapsladdergeos)
0207 {
0208 PHG4CylinderGeomContainer::ConstRange begin_end =
0209 mapsladdergeos->get_begin_end();
0210 PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
0211 for (; miter != begin_end.second; ++miter)
0212 {
0213 PHG4CylinderGeom *geo = miter->second;
0214
0215 //if(geo->get_layer() > (int) _radii.size() ) continue;
0216
0217 // if (Verbosity() >= 2)
0218 // geo->identify();
0219
0220 //TODO
0221 _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
0222 geo->get_radius();
0223
0224 /*if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
0225 {
0226 _radii[_layer_ilayer_map[geo->get_layer()]] = geo->get_radius();
0227 }*/
0228 }
0229 }
0230 return Fun4AllReturnCodes::EVENT_OK;
0231 }
0232
0233 int PHHybridSeeding::Process(PHCompositeNode */*topNode*/)
0234 {
0235 // wipe previous vertex coordinates
0236 _vertex_x.clear();
0237 _vertex_y.clear();
0238 _vertex_z.clear();
0239 _vertex_xerr.clear();
0240 _vertex_yerr.clear();
0241 _vertex_zerr.clear();
0242 _vertex_ids.clear();
0243 // fill new vertex coordinates
0244 for(map<unsigned int, SvtxVertex*>::iterator iter = _vertex_map->begin(); iter != _vertex_map->end(); ++iter)
0245 {
0246 SvtxVertex* v = dynamic_cast<SvtxVertex*>(iter->second->CloneMe());
0247 _vertex_x.push_back(v->get_x());
0248 _vertex_y.push_back(v->get_y());
0249 _vertex_z.push_back(v->get_z());
0250 _vertex_xerr.push_back(sqrt(v->get_error(0,0)));
0251 _vertex_yerr.push_back(sqrt(v->get_error(1,1)));
0252 _vertex_zerr.push_back(sqrt(v->get_error(2,2)));
0253 _vertex_ids.push_back(v->get_id());
0254 }
0255 if(Verbosity()>1) cout << "vertices:\n";
0256 for(size_t i=0;i<_vertex_x.size();i++)
0257 {
0258 if(Verbosity()>1) cout << "(" << _vertex_x[i] << "," << _vertex_y[i] << "," << _vertex_z[i] << ")\n";
0259 }
0260 vector<vector<double>> kdhits(PHTpcTrackerUtil::convert_clusters_to_hits(_cluster_map,_hitsets));
0261 vector<vector<double> > unused_hits;
0262 vector<vector<vector<double> > > kdtracks;
0263
0264 bool print_stats = (Verbosity()>0);
0265
0266 kdtracks = kdfinder::find_tracks_iterative<double>(kdhits, unused_hits,
0267 _search_radius1, /* max distance in cm*/
0268 _search_angle1, /* triplet angle */
0269 _min_track_size1, /* min hits to keep track */
0270 // first iteration
0271 _search_radius2,
0272 _search_angle2,
0273 _min_track_size2, // second iteration params
0274 _nthreads,
0275 print_stats);
0276 if(Verbosity()>0) cout << "n_kdtracks: " << kdtracks.size() << "\n";
0277 vector<keylist> clusterLists;
0278 for(auto track : kdtracks)
0279 {
0280 keylist k;
0281 for(auto kdhit : track)
0282 {
0283 // see PHTpcTracker/PHTpcTrackerUtil.cc; this recovers the cluster key, apparently
0284 k.push_back(*((int64_t*)&kdhit[3]));
0285 }
0286 clusterLists.push_back(k);
0287 }
0288 for(auto& clist : clusterLists)
0289 {
0290 if(Verbosity()>1) cout << "front layer: " << (int)TrkrDefs::getLayer(clist.front()) << " back layer: " << (int)TrkrDefs::getLayer(clist.back());
0291 if(clist.size()>1 && ((int)TrkrDefs::getLayer(clist.front()))<((int)TrkrDefs::getLayer(clist.back())))
0292 {
0293 if(Verbosity()>1) cout << "reversing\n";
0294 std::reverse(clist.begin(),clist.end());
0295 }
0296 if(Verbosity()>1) cout << "final layer order:\n";
0297 for(auto c : clist) if(Verbosity()>1) cout << (int)TrkrDefs::getLayer(c) << endl;
0298 }
0299 for(auto clist : clusterLists)
0300 {
0301 if(Verbosity()>1) cout << "layers:\n";
0302 for(auto c : clist)
0303 {
0304 if(Verbosity()>1) cout << (int)TrkrDefs::getLayer(c) << endl;
0305 }
0306 }
0307 vector<SvtxTrack_v3> seeds = fitter->ALICEKalmanFilter(clusterLists,false);
0308 if(Verbosity()>0) cout << "nseeds: " << seeds.size() << "\n";
0309 publishSeeds(seeds);
0310 return Fun4AllReturnCodes::EVENT_OK;
0311 }
0312
0313 void PHHybridSeeding::publishSeeds(vector<SvtxTrack_v3> seeds)
0314 {
0315 for(size_t i=0;i<seeds.size();i++)
0316 {
0317 _track_map->insert(&(seeds[i]));
0318 }
0319 }
0320
0321 int PHHybridSeeding::Setup(PHCompositeNode *topNode)
0322 {
0323 if(Verbosity()>0) cout << "Called Setup" << endl;
0324 if(Verbosity()>0) cout << "topNode:" << topNode << endl;
0325 PHTrackSeeding::Setup(topNode);
0326
0327 _vertex_map = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0328 if (!_vertex_map)
0329 {
0330 cerr << PHWHERE << " ERROR: Can't find SvtxVertexMap." << endl;
0331 return Fun4AllReturnCodes::ABORTEVENT;
0332 }
0333
0334 auto surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
0335 if(!surfmaps)
0336 {
0337 std::cout << "No Acts surface maps, bailing." << std::endl;
0338 return Fun4AllReturnCodes::ABORTEVENT;
0339 }
0340
0341 auto tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,"ActsTrackingGeometry");
0342 if(!tGeometry)
0343 {
0344 std::cout << "No Acts tracking geometry, exiting." << std::endl;
0345 return Fun4AllReturnCodes::ABORTEVENT;
0346 }
0347
0348 InitializeGeometry(topNode);
0349 fitter = std::make_shared<ALICEKF>(topNode,_cluster_map,surfmaps, tGeometry,
0350 _fieldDir,_min_fit_track_size,_max_sin_phi,Verbosity());
0351 return Fun4AllReturnCodes::EVENT_OK;
0352 }
0353
0354 int PHHybridSeeding::End()
0355 {
0356 if(Verbosity()>0) cout << "Called End " << endl;
0357 return Fun4AllReturnCodes::EVENT_OK;
0358 }
