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