Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHRaveVertexing.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHRaveVertexing.cc
1 
8 #include "PHRaveVertexing.h"
9 
12 #include <trackbase_historic/SvtxTrackState.h> // for SvtxTrackState
15 
17 
18 #include <phgenfit/Fitter.h>
19 
20 #include <phfield/PHFieldUtility.h>
21 
22 #include <phgeom/PHGeomUtility.h>
23 
25 #include <fun4all/SubsysReco.h> // for SubsysReco
26 
27 #include <phool/PHCompositeNode.h>
28 #include <phool/PHIODataNode.h>
29 #include <phool/PHNode.h> // for PHNode
30 #include <phool/PHNodeIterator.h>
31 #include <phool/PHObject.h> // for PHObject
32 #include <phool/PHTimer.h>
33 #include <phool/getClass.h>
34 #include <phool/phool.h>
35 
36 #include <GenFit/FitStatus.h> // for FitStatus
37 #include <GenFit/GFRaveTrackParameters.h> // for GFRaveTrackParameters
38 #include <GenFit/GFRaveVertex.h>
39 #include <GenFit/GFRaveVertexFactory.h>
40 #include <GenFit/KalmanFittedStateOnPlane.h> // for KalmanFittedStateOn...
41 #include <GenFit/KalmanFitterInfo.h>
42 #include <GenFit/MeasuredStateOnPlane.h>
43 #include <GenFit/RKTrackRep.h>
44 #include <GenFit/Track.h>
45 #include <GenFit/TrackPoint.h> // for TrackPoint
46 
47 #include <TMatrixDSymfwd.h> // for TMatrixDSym
48 #include <TMatrixTSym.h> // for TMatrixTSym
49 #include <TMatrixTUtils.h> // for TMatrixTRow
50 #include <TVector3.h>
51 
52 #include <iostream>
53 #include <map>
54 #include <memory>
55 #include <utility>
56 #include <vector>
57 
58 class PHField;
59 class TGeoManager;
60 namespace genfit { class AbsTrackRep; }
61 
62 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
63 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
64 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
65 
66 //#define _DEBUG_
67 
68 using namespace std;
69 
70 /*
71  * Constructor
72  */
74  : SubsysReco(name)
75  , _over_write_svtxvertexmap(false)
76  , _svtxvertexmaprefit_node_name("SvtxVertexMapRefit")
77  , _fitter(nullptr)
78  , _primary_pid_guess(211)
79  , _vertex_min_ndf(20)
80  , _vertex_finder(nullptr)
81  , _vertexing_method("avf-smoothing:1")
82  , _trackmap(nullptr)
83  , _vertexmap(nullptr)
84  , _vertexmap_refit(nullptr)
85  , _t_translate(nullptr)
86  , _t_rave(nullptr)
87 {
88  Verbosity(0);
89 
90  _event = 0;
91 }
92 
93 /*
94  * Init
95  */
97 {
99 }
100 
101 /*
102  * Init run
103  */
105 {
106  CreateNodes(topNode);
107 
108  TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
109  PHField* field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
110 
111  //_fitter = new PHGenFit::Fitter("sPHENIX_Geo.root","sPHENIX.2d.root", 1.4 / 1.5);
112  _fitter = PHGenFit::Fitter::getInstance(tgeo_manager,
113  field, "DafRef",
114  "RKTrackRep", false);
115  if (!_fitter)
116  {
117  std::cout << PHWHERE << " PHGenFit::Fitter::getInstance returned nullptr"
118  << std::endl;
120  }
122 
123  if (!_fitter)
124  {
125  cerr << PHWHERE << endl;
127  }
128 
130  if (!_vertex_finder)
131  {
132  std::cout << PHWHERE << " genfit::GFRaveVertexFactory returned null ptr" << std::endl;
134  }
136  //_vertex_finder->setBeamspot();
137 
138 
139  _t_translate = new PHTimer("_t_translate");
140  _t_translate->stop();
141 
142  _t_rave = new PHTimer("_t_rave");
143  _t_rave->stop();
144 
146 }
154 {
155  _event++;
156 
157  if (Verbosity() > 1)
158  std::cout << PHWHERE << "Events processed: " << _event << std::endl;
159 
160  GetNodes(topNode);
161 
163  GenFitTrackMap gf_track_map;
164  vector<genfit::Track*> gf_tracks;
165  if (Verbosity() > 1) _t_translate->restart();
166  for (SvtxTrackMap::Iter iter = _trackmap->begin(); iter != _trackmap->end();
167  ++iter)
168  {
169  SvtxTrack* svtx_track = iter->second;
170  if (!svtx_track)
171  continue;
172 
173  if (!(svtx_track->get_ndf() >= _vertex_min_ndf))
174  continue;
175 
176  // require MVTX association
177  if(_nmvtx_required > 0)
178  {
179  unsigned int nmvtx = 0;
180  for(auto clusit = svtx_track->begin_cluster_keys(); clusit != svtx_track->end_cluster_keys(); ++clusit)
181  {
182  if(TrkrDefs::getTrkrId(*clusit) == TrkrDefs::mvtxId )
183  {
184  nmvtx++;
185  }
186  if(nmvtx >= _nmvtx_required) break;
187  }
188  if(nmvtx < _nmvtx_required) continue;
189  if(Verbosity() > 1) std::cout << " track " << iter->first << " has nmvtx at least " << nmvtx << std::endl;
190  }
191 
192  //auto genfit_track = shared_ptr<genfit::Track> (TranslateSvtxToGenFitTrack(svtx_track));
193  auto genfit_track = TranslateSvtxToGenFitTrack(svtx_track);
194  if (!genfit_track)
195  continue;
196  gf_track_map.insert({genfit_track, iter->first});
197  gf_tracks.push_back(const_cast<genfit::Track*>(genfit_track));
198  }
199  if (Verbosity() > 1) _t_translate->stop();
200 
201  if (Verbosity() > 1) _t_rave->restart();
202  vector<genfit::GFRaveVertex*> rave_vertices;
203  if (gf_tracks.size() >= 2)
204  {
205  try
206  {
207  _vertex_finder->findVertices(&rave_vertices, gf_tracks);
208  }
209  catch (...)
210  {
211  if (Verbosity() > 1)
212  std::cout << PHWHERE << "GFRaveVertexFactory::findVertices failed!";
213  }
214  }
215  if (Verbosity() > 1) _t_rave->stop();
216  FillSvtxVertexMap(rave_vertices, gf_track_map);
217 
218  for (auto iter : gf_track_map) delete iter.first;
219 
220  if (Verbosity() > 1)
221  {
222  _vertexmap->identify();
223 
224  std::cout << "=============== Timers: ===============" << std::endl;
225  std::cout << "Event: " << _event << std::endl;
226  std::cout << "Translate: " << _t_translate->get_accumulated_time() / 1000. << " sec" << std::endl;
227  std::cout << "RAVE: " << _t_rave->get_accumulated_time() / 1000. << " sec" << std::endl;
228  std::cout << "=======================================" << std::endl;
229  }
230 
232 }
233 
234 /*
235  * End
236  */
238 {
240 }
241 
242 /*
243  * dtor
244  */
246 {
247  delete _fitter;
248  delete _vertex_finder;
249 }
250 
252 {
253  // create nodes...
254  PHNodeIterator iter(topNode);
255 
256  PHCompositeNode* dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
257  "PHCompositeNode", "DST"));
258  if (!dstNode)
259  {
260  cerr << PHWHERE << "DST Node missing, doing nothing." << endl;
262  }
263  PHNodeIterator iter_dst(dstNode);
264 
265  // Create the SVTX node
266  PHCompositeNode* tb_node = dynamic_cast<PHCompositeNode*>(iter_dst.findFirst(
267  "PHCompositeNode", "SVTX"));
268  if (!tb_node)
269  {
270  tb_node = new PHCompositeNode("SVTX");
271  dstNode->addNode(tb_node);
272  if (Verbosity() > 0)
273  cout << "SVTX node added" << endl;
274  }
275 
277  {
279  PHIODataNode<PHObject>* vertexes_node = new PHIODataNode<PHObject>(
280  _vertexmap_refit, _svtxvertexmaprefit_node_name.c_str(), "PHObject");
281  tb_node->addNode(vertexes_node);
282  if (Verbosity() > 0)
283  cout << "Svtx/SvtxVertexMapRefit node added" << endl;
284  }
285  else if (!findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap"))
286  {
288  PHIODataNode<PHObject>* vertexes_node = new PHIODataNode<PHObject>(
289  _vertexmap, "SvtxVertexMap", "PHObject");
290  tb_node->addNode(vertexes_node);
291  if (Verbosity() > 0)
292  cout << "Svtx/SvtxVertexMap node added" << endl;
293  }
294 
296 }
297 
298 /*
299  * GetNodes():
300  * Get all the all the required nodes off the node tree
301  */
303 {
304  //DST objects
305  // Input Svtx Tracks
306  _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
307  if (!_trackmap && _event < 2)
308  {
309  cout << PHWHERE << " SvtxTrackMap node not found on node tree"
310  << endl;
312  }
313 
314  // Input Svtx Vertices
315  _vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
316  if (!_vertexmap && _event < 2)
317  {
318  cout << PHWHERE << " SvtxVertexrMap node not found on node tree"
319  << endl;
321  }
322 
323  // Output Svtx Vertices
325  {
326  _vertexmap_refit = findNode::getClass<SvtxVertexMap>(topNode,
328  if (!_vertexmap_refit && _event < 2)
329  {
330  cout << PHWHERE << " SvtxVertexMapRefit node not found on node tree"
331  << endl;
333  }
334  }
335 
337 }
338 
339 /*
340  * Fill SvtxVertexMap from GFRaveVertexes and Tracks
341  */
343  const std::vector<genfit::GFRaveVertex*>& rave_vertices,
344  const GenFitTrackMap& gf_track_map)
345 {
347  {
348  _vertexmap->clear();
349  }
350 
351  // for (unsigned int ivtx = 0; ivtx < rave_vertices.size(); ++ivtx) {
352  // genfit::GFRaveVertex* rave_vtx = rave_vertices[ivtx];
353 
354  for (genfit::GFRaveVertex* rave_vtx : rave_vertices)
355  {
356  if (!rave_vtx)
357  {
358  cerr << PHWHERE << endl;
359  return false;
360  }
361 
362  std::shared_ptr<SvtxVertex> svtx_vtx(new SvtxVertex_v1());
363 
364  svtx_vtx->set_chisq(rave_vtx->getChi2());
365  svtx_vtx->set_ndof(rave_vtx->getNdf());
366  svtx_vtx->set_position(0, rave_vtx->getPos().X());
367  svtx_vtx->set_position(1, rave_vtx->getPos().Y());
368  svtx_vtx->set_position(2, rave_vtx->getPos().Z());
369 
370  for (int i = 0; i < 3; i++)
371  for (int j = 0; j < 3; j++)
372  svtx_vtx->set_error(i, j, rave_vtx->getCov()[i][j]);
373 
374  for (unsigned int i = 0; i < rave_vtx->getNTracks(); i++)
375  {
376  //TODO improve speed
377  const genfit::Track* rave_track =
378  rave_vtx->getParameters(i)->getTrack();
379  // for(auto iter : gf_track_map) {
380  // if (iter.second == rave_track)
381  // svtx_vtx->insert_track(iter.first);
382  // }
383  auto iter = gf_track_map.find(rave_track);
384  if (iter != gf_track_map.end())
385  svtx_vtx->insert_track(iter->second);
386  }
387 
389  {
390  if (_vertexmap)
391  {
392  _vertexmap->insert_clone(svtx_vtx.get());
393  }
394  else
395  {
396  LogError("!_vertexmap");
397  }
398  }
399  else
400  {
401  if (_vertexmap_refit)
402  {
403  _vertexmap_refit->insert_clone(svtx_vtx.get());
404  }
405  else
406  {
407  LogError("!_vertexmap_refit");
408  }
409  }
410 
411 #ifdef _DEBUG_
412  cout << __LINE__ << endl;
413  svtx_vtx->identify();
414 #endif
415 
416  } //loop over RAVE vertices
417 
418  return true;
419 }
420 
422 {
423  try
424  {
425  // The first state is extracted to PCA, second one is the one with measurement
426  SvtxTrackState* svtx_state(nullptr);
427  //SvtxTrackState* svtx_state = (svtx_track->begin_states())->second;
428 
429  if (svtx_track->begin_states() == svtx_track->end_states())
430  {
431  LogDebug("TranslateSvtxToGenFitTrack no state in track!");
432  return nullptr;
433  }
434  else if (++(svtx_track->begin_states()) == svtx_track->end_states())
435  {
436  // only one state in track
437  svtx_state = (svtx_track->begin_states())->second;
438  }
439  else
440  {
441  // multiple state in track
442  // The first state is extracted to PCA, second one is the one with measurement
443  svtx_state = (++(svtx_track->begin_states()))->second;
444  }
445 
446  if (!svtx_state)
447  {
448  LogDebug("TranslateSvtxToGenFitTrack invalid state found on track!");
449  return nullptr;
450  }
451 
452  TVector3 pos(svtx_state->get_x(), svtx_state->get_y(), svtx_state->get_z());
453  TVector3 mom(svtx_state->get_px(), svtx_state->get_py(), svtx_state->get_pz());
454  TMatrixDSym cov(6);
455  for (int i = 0; i < 6; ++i)
456  {
457  for (int j = 0; j < 6; ++j)
458  {
459  cov[i][j] = svtx_state->get_error(i, j);
460  }
461  }
462 
463 #ifdef _DEBUG_
464  {
465  cout << "DEBUG" << __LINE__ << endl;
466  cout << "path length: " << svtx_state->get_pathlength() << endl;
467  cout << "radius: " << pos.Perp() << endl;
468  cout << "DEBUG: " << __LINE__ << endl;
469  for (int i = 0; i < 6; ++i)
470  {
471  for (int j = 0; j < 6; ++j)
472  {
473  cout << svtx_state->get_error(i, j) << "\t";
474  }
475  cout << endl;
476  }
477 
478  cov.Print();
479  }
480 
481 #endif
482 
484  genfit::Track* genfit_track = new genfit::Track(rep, TVector3(0, 0, 0), TVector3(0, 0, 0));
485 
487  fs->setCharge(svtx_track->get_charge());
488  fs->setChi2(svtx_track->get_chisq());
489  fs->setNdf(svtx_track->get_ndf());
490  fs->setIsFitted(true);
491  fs->setIsFitConvergedFully(true);
492 
493  genfit_track->setFitStatus(fs, rep);
494 
495  genfit::TrackPoint* tp = new genfit::TrackPoint(genfit_track);
496 
498  tp->setFitterInfo(fi);
499 
501  ms->setPosMomCov(pos, mom, cov);
502 #ifdef _DEBUG_
503  {
504  cout << "DEBUG: " << __LINE__ << endl;
505  ms->Print();
506  cout << "Orig: " << __LINE__ << endl;
507  cov.Print();
508  cout << "Translate: " << __LINE__ << endl;
509  ms->get6DCov().Print();
510  }
511 #endif
513 
514  //< Acording to the special order of using the stored states
515  fi->setForwardUpdate(kfs);
516 
517  genfit_track->insertPoint(tp);
518 
519 #ifdef _DEBUG_
520 // {
521 // cout << "DEBUG" << __LINE__ << endl;
522 // TVector3 pos, mom;
523 // TMatrixDSym cov;
524 // genfit_track->getFittedState().getPosMomCov(pos, mom, cov);
525 // pos.Print();
526 // mom.Print();
527 // cov.Print();
528 // }
529 #endif
530 
531  return genfit_track;
532  }
533  catch (...)
534  {
535  LogDebug("TranslateSvtxToGenFitTrack failed!");
536  }
537 
538  return nullptr;
539 }