ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TrackFastSim.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TrackFastSim.cc
1 
8 #include "PHG4TrackFastSim.h"
9 
10 #include <phgenfit/Fitter.h>
11 #include <phgenfit/Measurement.h> // for Measurement
12 #include <phgenfit/PlanarMeasurement.h>
13 #include <phgenfit/SpacepointMeasurement.h>
14 #include <phgenfit/Track.h>
15 #include <phgeom/PHGeomUtility.h>
16 
23 #include <trackbase_historic/SvtxVertex.h> // for SvtxVertex
24 #include <trackbase_historic/SvtxVertexMap.h> // for SvtxVertexMap
27 
28 #include <calobase/RawTowerGeom.h>
29 #include <calobase/RawTowerGeomContainer.h>
30 
31 #include <phparameter/PHParameters.h>
32 
33 #include <phfield/PHFieldUtility.h>
34 
35 #include <g4main/PHG4Hit.h>
36 #include <g4main/PHG4HitContainer.h> // for PHG4HitContainer
37 #include <g4main/PHG4Particle.h>
39 #include <g4main/PHG4VtxPoint.h>
40 
42 #include <fun4all/SubsysReco.h> // for SubsysReco
43 
44 #include <phool/PHCompositeNode.h>
45 #include <phool/PHIODataNode.h>
46 #include <phool/PHNode.h> // for PHNode
47 #include <phool/PHNodeIterator.h>
48 #include <phool/PHObject.h> // for PHObject
49 #include <phool/PHRandomSeed.h>
50 #include <phool/getClass.h>
51 #include <phool/phool.h>
52 
53 #include <GenFit/AbsMeasurement.h>
54 #include <GenFit/EventDisplay.h>
55 #include <GenFit/MeasuredStateOnPlane.h>
56 #include <GenFit/RKTrackRep.h>
57 
58 #include <GenFit/FitStatus.h> // for FitStatus
59 #include <GenFit/GFRaveTrackParameters.h> // for GFRaveTrackParameters
60 #include <GenFit/GFRaveVertex.h>
61 #include <GenFit/GFRaveVertexFactory.h>
62 #include <GenFit/Track.h>
63 
64 #include <TMatrixDSymfwd.h> // for TMatrixDSym
65 #include <TMatrixTSym.h> // for TMatrixTSym
66 #include <TMatrixTUtils.h> // for TMatrixTRow
67 #include <TSystem.h>
68 #include <TVector3.h> // for TVector3, operator*
69 #include <TVectorDfwd.h> // for TVectorD
70 #include <TVectorT.h> // for TVectorT
71 
72 #include <gsl/gsl_randist.h>
73 #include <gsl/gsl_rng.h>
74 
75 #include <cassert> // for assert
76 #include <cmath>
77 #include <iostream> // for operator<<, basic_...
78 #include <map>
79 #include <memory> // for unique_ptr, alloca...
80 #include <utility>
81 
82 class PHField;
83 class TGeoManager;
84 namespace genfit
85 {
86  class AbsTrackRep;
87 } // namespace genfit
88 
89 #define LogDebug(exp) \
90  if (Verbosity()) std::cout << "PHG4TrackFastSim (DEBUG): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
91 #define LogError(exp) \
92  std::cout << "PHG4TrackFastSim (ERROR): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
93 #define LogWarning(exp) \
94  std::cout << "PHG4TrackFastSim (WARNING): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
95 
96 using namespace std;
97 
98 // names of our implemented calorimeters where the projections are done
99 // at 1/2 of their depth, not at the surface
100 // this is used to avoid user added projections with identical names
101 set<string> reserved_cylinder_projection_names{"CEMC", "HCALIN", "HCALOUT", "BECAL"};
102 set<string> reserved_zplane_projection_names{"FEMC", "FHCAL", "EEMC", "EHCAL", "LFHCAL"};
103 
105  : SubsysReco(name)
106  , m_Fitter(nullptr)
107  , m_RaveVertexFactory(nullptr)
108  , m_TruthContainer(nullptr)
109  , m_SvtxTrackMapOut(nullptr)
110  , m_SvtxVertexMap(nullptr)
111  , m_SubTopnodeName("SVTX")
112  , m_TrackmapOutNodeName("SvtxTrackMap")
113  , m_VertexingMethod("kalman-smoothing:1")
114  , m_FitAlgoName("DafRef") // was ("KalmanFitterRefTrack")
115  , m_VertexMinNdf(10.)
116  , m_VertexXYResolution(50E-4)
117  , m_VertexZResolution(50E-4)
118  , m_EventCnt(-1)
119  , m_PrimaryAssumptionPid(211)
120  , m_SmearingFlag(true)
121  , m_DoEvtDisplayFlag(false)
122  , m_UseVertexInFittingFlag(true)
123  , m_PrimaryTrackingFlag(1)
124  , m_DoVertexingFlag(false)
125 {
126  unsigned int seed = PHRandomSeed(); // fixed seed is handled in this funtcion
127  m_RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
128  gsl_rng_set(m_RandomGenerator, seed);
129 
130  m_Parameter = new PHParameters(Name());
131 }
132 
134 {
135  delete m_Fitter;
136  delete m_RaveVertexFactory;
137  gsl_rng_free(m_RandomGenerator);
138 }
139 
141 {
142  m_EventCnt = -1;
143 
144  int ret = CreateNodes(topNode);
145  if (ret != Fun4AllReturnCodes::EVENT_OK)
146  {
147  return ret;
148  }
149  ret = GetNodes(topNode);
150  if (ret != Fun4AllReturnCodes::EVENT_OK)
151  {
152  return ret;
153  }
154  //
155  //
156  TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
157  PHField* field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
158 
160  field, m_FitAlgoName, "RKTrackRep",
162 
163  if (!m_Fitter)
164  {
165  cerr << PHWHERE << endl;
167  }
168 
170 
171  // tower geometry for track states
172 
173  for (map<string, pair<int, double>>::iterator iter = m_ProjectionsMap.begin(); iter != m_ProjectionsMap.end(); ++iter)
174  {
175  if (isfinite(iter->second.second))
176  {
177  continue;
178  }
179  switch (iter->second.first)
180  {
181  case DETECTOR_TYPE::Cylinder:
182  {
183  string nodename = "TOWERGEOM_" + iter->first;
184  RawTowerGeomContainer* geo = findNode::getClass<RawTowerGeomContainer>(topNode, nodename);
185  if (geo)
186  {
187  iter->second.second = geo->get_radius();
188  }
189  break;
190  }
191  case DETECTOR_TYPE::Vertical_Plane:
192  {
193  string towergeonodename = "TOWERGEOM_" + iter->first;
194  RawTowerGeomContainer* towergeo = findNode::getClass<RawTowerGeomContainer>(topNode, towergeonodename);
195  if (!towergeo)
196  {
197  cout << PHWHERE << " ERROR: Can't find node " << towergeonodename << endl;
199  }
200 
201  // Grab the first tower, us it to get the
202  // location along the beamline
204  RawTowerGeomContainer::ConstIterator twr_iter = twr_range.first;
205  RawTowerGeom* temp_geo = twr_iter->second;
206 
207  //Changed by Barak on 12/10/19
208  iter->second.second = temp_geo->get_center_z();
209  break;
210  }
211  default:
212  cout << "invalid state reference: " << iter->second.first << endl;
213  gSystem->Exit(1);
214  }
215  }
216 
217  if (m_DoVertexingFlag)
218  {
219  m_RaveVertexFactory = new genfit::GFRaveVertexFactory(Verbosity(), true);
220  //m_RaveVertexFactory->setMethod("kalman-smoothing:1"); //! kalman-smoothing:1 is the defaul method
222  //m_RaveVertexFactory->setBeamspot();
223 
224  //m_RaveVertexFactory = new PHRaveVertexFactory(Verbosity());
225 
226  if (!m_RaveVertexFactory)
227  {
228  cout << PHWHERE << " no Vertex Finder" << endl;
230  }
231  }
233 }
234 
236 {
238  {
240  }
241 
243 }
244 
246 {
247  m_EventCnt++;
248 
249  if (Verbosity() >= 2)
250  std::cout << "PHG4TrackFastSim::process_event: " << m_EventCnt << ".\n";
251 
252  // if(_clustermap_out)
253  // _clustermap_out->empty();
254  // else {
255  // LogError("_clustermap_out not found!");
256  // return Fun4AllReturnCodes::ABORTRUN;
257  // }
258 
259  if (not m_SvtxTrackMapOut)
260  {
261  LogError("m_SvtxTrackMapOut not found!");
263  }
264 
265  vector<PHGenFit::Track*> rf_tracks;
266 
268  TVector3 vtxPoint(truthVtx->get_x(), truthVtx->get_y(), truthVtx->get_z());
269  // Smear the vertex ONCE for all particles in the event
270  vtxPoint.SetX(vtxPoint.x() + gsl_ran_gaussian(m_RandomGenerator, m_VertexXYResolution));
271  vtxPoint.SetY(vtxPoint.y() + gsl_ran_gaussian(m_RandomGenerator, m_VertexXYResolution));
272  vtxPoint.SetZ(vtxPoint.z() + gsl_ran_gaussian(m_RandomGenerator, m_VertexZResolution));
273 
276  {
277  // Tracking for primaries only
279  }
280  else
281  {
282  // Check ALL particles
283  itr_range = m_TruthContainer->GetParticleRange();
284  }
285 
286  GenFitTrackMap gf_track_map;
287  // Now we can loop over the particles
288 
289  for (PHG4TruthInfoContainer::ConstIterator itr = itr_range.first;
290  itr != itr_range.second; ++itr)
291  {
292  PHG4Particle* particle = itr->second;
293 
294  TVector3 seed_pos(vtxPoint.x(), vtxPoint.y(), vtxPoint.z());
295  TVector3 seed_mom(0, 0, 0);
296  TMatrixDSym seed_cov(6);
297 
299  std::vector<PHGenFit::Measurement*> measurements;
300 
301  PHGenFit::Measurement* vtx_meas = nullptr;
302 
304  {
305  vtx_meas = VertexMeasurement(TVector3(vtxPoint.x(),
306  vtxPoint.y(),
307  vtxPoint.z()),
310  measurements.push_back(vtx_meas);
311  }
312 
313  unique_ptr<SvtxTrack> svtx_track_out(new SvtxTrack_FastSim_v1());
314  PseudoPatternRecognition(particle,
315  measurements,
316  svtx_track_out.get(),
317  seed_pos, seed_mom,
318  seed_cov);
319 
320  if (measurements.size() < 3)
321  {
322  if (Verbosity() >= 2)
323  {
324  //LogWarning("measurements.size() < 3");
325  std::cout << "event: " << m_EventCnt << " : measurements.size() < 3"
326  << "\n";
327  }
328  // Delete the measurements
329  // We need to also delete the underlying genfit::AbsMeasurement object
330  for (unsigned int im = 0; im < measurements.size(); im++)
331  {
332  delete measurements[im]->getMeasurement();
333  delete measurements[im];
334  }
335  continue;
336  }
337 
339 
347  //int pid = 13; //
348  //SMART(genfit::AbsTrackRep) rep = NEW(genfit::RKTrackRep)(pid);
349  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(m_PrimaryAssumptionPid);
350 
351  //rep->setDebugLvl(1); //DEBUG
352 
354 
355  PHGenFit::Track* track = new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov);
356 
357  // NOTE: We need to keep a list of tracks so they can
358  // all be cleanly deleted at the end
359  rf_tracks.push_back(track);
360 
361  //LogDEBUG;
363  track->addMeasurements(measurements);
364 
365  //LogDEBUG;
367  int fitting_err = m_Fitter->processTrack(track, false);
368 
369  if (fitting_err != 0)
370  {
371  if (Verbosity() >= 2)
372  {
373  //LogWarning("measurements.size() < 3");
374  std::cout << "event: " << m_EventCnt
375  << " : fitting_err != 0, next track."
376  << "\n";
377  }
378  continue;
379  }
380 
381  TVector3 vtx(vtxPoint.x(), vtxPoint.y(), vtxPoint.z());
382  bool track_made = MakeSvtxTrack(svtx_track_out.get(), track,
383  particle->get_track_id(),
384  measurements.size(), vtx);
385  if (Verbosity() > 1)
386  {
387  svtx_track_out->identify();
388  }
389 
390  if (track_made)
391  {
392  // track -> output container
393 
394  const unsigned int track_id = m_SvtxTrackMapOut->insert(svtx_track_out.get())->get_id();
395  gf_track_map.insert({track->getGenFitTrack(), track_id});
396  }
397 
398  } // Loop all primary particles
399 
400  //vertex finding
401  if (m_DoVertexingFlag)
402  {
403  if (!m_RaveVertexFactory)
404  {
405  cout << __PRETTY_FUNCTION__ << "Failed to init vertex finder" << endl;
407  }
408  if (!m_SvtxVertexMap)
409  {
410  cout << __PRETTY_FUNCTION__ << "Failed to init vertex map" << endl;
412  }
413 
414  // genfit::GFRaveVertexFactory* m_RaveVertexFactory = new genfit::GFRaveVertexFactory(10, true);
415  // m_RaveVertexFactory->setMethod("kalman-smoothing:1");
416  // m_RaveVertexFactory->setBeamspot();
417 
418  vector<genfit::GFRaveVertex*> rave_vertices;
419  if (rf_tracks.size() >= 2)
420  {
421  try
422  {
423  vector<genfit::Track*> rf_gf_tracks;
424  for (std::vector<PHGenFit::Track*>::iterator it = rf_tracks.begin(); it != rf_tracks.end(); ++it)
425  {
426  genfit::Track* track = (*it)->getGenFitTrack();
427 
428  if (Verbosity())
429  {
430  TVector3 pos, mom;
431  TMatrixDSym cov;
432 
433  track->getFittedState().getPosMomCov(pos, mom, cov);
434 
435  cout << "Track getCharge = " << track->getFitStatus()->getCharge() << " getChi2 = " << track->getFitStatus()->getChi2() << " getNdf = " << track->getFitStatus()->getNdf() << endl;
436  pos.Print();
437  mom.Print();
438  cov.Print();
439  }
440  if (track->getFitStatus()->getNdf() > m_VertexMinNdf)
441  rf_gf_tracks.push_back(track);
442  }
443  m_RaveVertexFactory->findVertices(&rave_vertices, rf_gf_tracks);
444  }
445  catch (...)
446  {
447  if (Verbosity() > 1)
448  std::cout << PHWHERE << "GFRaveVertexFactory::findVertices failed!";
449  }
450  }
451 
452  if (Verbosity())
453  {
454  cout << __PRETTY_FUNCTION__ << __LINE__ << " rf_tracks = " << rf_tracks.size() << " rave_vertices = " << rave_vertices.size() << endl;
455  }
456  FillSvtxVertexMap(rave_vertices, gf_track_map);
457  }
458 
460  if (m_DoEvtDisplayFlag)
461  {
462  vector<genfit::Track*> rf_gf_tracks;
463  for (std::vector<PHGenFit::Track*>::iterator it = rf_tracks.begin(); it != rf_tracks.end(); ++it)
464  {
465  rf_gf_tracks.push_back((*it)->getGenFitTrack());
466  }
467  m_Fitter->getEventDisplay()->addEvent(rf_gf_tracks);
468  }
469  else
470  {
471  for (std::vector<PHGenFit::Track*>::iterator it = rf_tracks.begin(); it != rf_tracks.end(); ++it)
472  {
473  delete (*it);
474  }
475  rf_tracks.clear();
476  }
477 
478  // if(m_SvtxTrackMapOut->get(0)) {
479  // m_SvtxTrackMapOut->get(0)->identify();
480  // std::cout<<"DEBUG : "<< m_SvtxTrackMapOut->get(0)->get_px() <<"\n";
481  // std::cout<<"DEBUG : "<< m_SvtxTrackMapOut->get(0)->get_truth_track_id() <<"\n";
482  // }
483 
485 }
486 
487 /*
488  * Fill SvtxVertexMap from GFRaveVertexes and Tracks
489  */
491  const std::vector<genfit::GFRaveVertex*>& rave_vertices,
492  const GenFitTrackMap& gf_track_map)
493 {
494  for (genfit::GFRaveVertex* rave_vtx : rave_vertices)
495  {
496  if (!rave_vtx)
497  {
498  cerr << PHWHERE << endl;
499  return false;
500  }
501 
502  std::shared_ptr<SvtxVertex> svtx_vtx(new SvtxVertex_v1());
503 
504  svtx_vtx->set_chisq(rave_vtx->getChi2());
505  svtx_vtx->set_ndof(rave_vtx->getNdf());
506  svtx_vtx->set_position(0, rave_vtx->getPos().X());
507  svtx_vtx->set_position(1, rave_vtx->getPos().Y());
508  svtx_vtx->set_position(2, rave_vtx->getPos().Z());
509 
510  for (int i = 0; i < 3; i++)
511  {
512  for (int j = 0; j < 3; j++)
513  {
514  svtx_vtx->set_error(i, j, rave_vtx->getCov()[i][j]);
515  }
516  }
517 
518  for (unsigned int i = 0; i < rave_vtx->getNTracks(); i++)
519  {
520  //TODO improve speed
521  const genfit::Track* rave_track =
522  rave_vtx->getParameters(i)->getTrack();
523  // for(auto iter : gf_track_map) {
524  // if (iter.second == rave_track)
525  // svtx_vtx->insert_track(iter.first);
526  // }
527  auto iter = gf_track_map.find(rave_track);
528  if (iter != gf_track_map.end())
529  {
530  svtx_vtx->insert_track(iter->second);
531  }
532  }
533 
534  if (m_SvtxVertexMap)
535  {
536  m_SvtxVertexMap->insert_clone(svtx_vtx.get());
537  }
538  else
539  {
540  LogError("!m_SvtxVertexMap");
541  }
542 
543  } //loop over RAVE vertices
544 
545  return true;
546 }
547 
549 {
550  // create nodes...
551  PHNodeIterator iter(topNode);
552 
553  PHCompositeNode* dstNode = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
554  if (!dstNode)
555  {
556  cout << PHWHERE << " DST Node missing, doing nothing." << endl;
558  }
559  PHNodeIterator iter_dst(dstNode);
560 
561  // Create the FGEM node
562  PHCompositeNode* tb_node = dynamic_cast<PHCompositeNode*>(iter_dst.findFirst(
563  "PHCompositeNode", m_SubTopnodeName));
564  if (!tb_node)
565  {
566  tb_node = new PHCompositeNode(m_SubTopnodeName);
567  dstNode->addNode(tb_node);
568  if (Verbosity() > 0)
569  {
570  cout << m_SubTopnodeName << " node added" << endl;
571  }
572  }
573 
574  // _clustermap_out = new SvtxClusterMap_v1;
575  //
576  // PHIODataNode<PHObject>* clusters_node = new PHIODataNode<PHObject>(
577  // _clustermap_out, _clustermap_out_name, "PHObject");
578  // tb_node->addNode(clusters_node);
579  // if (Verbosity() > 0)
580  // cout << _clustermap_out_name <<" node added" << endl;
581 
582  m_SvtxTrackMapOut = findNode::getClass<SvtxTrackMap>(topNode, m_TrackmapOutNodeName);
583  if (!m_SvtxTrackMapOut)
584  {
586 
588  tb_node->addNode(tracks_node);
589  if (Verbosity() > 0)
590  {
591  cout << m_TrackmapOutNodeName << " node added" << endl;
592  }
593  }
594 
595  m_SvtxVertexMap = findNode::getClass<SvtxVertexMap_v1>(topNode, "SvtxVertexMap");
596  if (not m_SvtxVertexMap)
597  {
598  m_SvtxVertexMap = new SvtxVertexMap_v1;
599  PHIODataNode<PHObject>* vertexes_node = new PHIODataNode<PHObject>(m_SvtxVertexMap, "SvtxVertexMap", "PHObject");
600  tb_node->addNode(vertexes_node);
601  if (Verbosity() > 0)
602  {
603  cout << "Svtx/SvtxVertexMap node added" << endl;
604  }
605  }
606 
608 }
609 
611 {
612  assert(m_Parameter);
613  PHNodeIterator iter(topNode);
614 
615  //DST objects
616  //Truth container
617  m_TruthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
618  if (!m_TruthContainer)
619  {
620  cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree"
621  << endl;
623  }
624 
625  for (unsigned int i = 0; i < m_PHG4HitsNames.size(); i++)
626  {
627  PHG4HitContainer* phg4hit = findNode::getClass<PHG4HitContainer>(topNode, m_PHG4HitsNames[i]);
628  if (!phg4hit)
629  {
630  cout << PHWHERE << m_PHG4HitsNames[i]
631  << " node not found on node tree" << endl;
633  }
634 
635  if (Verbosity() > 0)
636  {
637  cout << "PHG4TrackFastSim::GetNodes - node added: " << m_PHG4HitsNames[i] << endl;
638  }
639  m_PHG4HitContainer.push_back(phg4hit);
640 
642  }
643 
644  // Update Kalman filter parameter node
645  m_Parameter->set_string_param("SubTopnodeName", m_SubTopnodeName);
646  m_Parameter->set_string_param("TrackmapOutNodeName", m_TrackmapOutNodeName);
647  m_Parameter->set_string_param("VertexingMethod", m_VertexingMethod);
648  m_Parameter->set_string_param("FitAlgoName", m_FitAlgoName);
649 
650  PHCompositeNode* runNode = static_cast<PHCompositeNode*>(iter.findFirst(
651  "PHCompositeNode", "RUN"));
652  if (!runNode)
653  {
654  std::cout << Name() << "::"
655  << "::" << __PRETTY_FUNCTION__
656  << "Run Node missing!" << std::endl;
657  throw std::runtime_error(
658  "Failed to find Run node in PHG4TrackFastSim::CreateNodes");
659  }
660  if (Verbosity())
661  {
662  cout << __PRETTY_FUNCTION__ << " : ";
664  }
665  m_Parameter->SaveToNodeTree(runNode, Name() + string("_Parameter"));
666 
667  //checks
668  assert(m_PHG4HitsNames.size() == m_PHG4HitContainer.size());
669  assert(m_phg4_detector_type.size() == m_PHG4HitContainer.size());
670  assert(m_phg4_detector_radres.size() == m_PHG4HitContainer.size());
671  assert(m_phg4_detector_phires.size() == m_PHG4HitContainer.size());
672  assert(m_phg4_detector_lonres.size() == m_PHG4HitContainer.size());
673  assert(m_phg4_detector_hitfindeff.size() == m_PHG4HitContainer.size());
674  assert(m_phg4_detector_noise.size() == m_PHG4HitContainer.size());
675 
676  // _clustermap_out = findNode::getClass<SvtxClusterMap>(topNode,
677  // _clustermap_out_name);
678  // if (!_clustermap_out && m_EventCnt < 2) {
679  // cout << PHWHERE << _clustermap_out_name << " node not found on node tree"
680  // << endl;
681  // return Fun4AllReturnCodes::ABORTEVENT;
682  // }
683 
684  m_SvtxTrackMapOut = findNode::getClass<SvtxTrackMap>(topNode, m_TrackmapOutNodeName);
685  if (!m_SvtxTrackMapOut && m_EventCnt < 2)
686  {
687  cout << PHWHERE << m_TrackmapOutNodeName
688  << " node not found on node tree" << endl;
690  }
691 
693 }
694 
696  std::vector<PHGenFit::Measurement*>& meas_out,
697  SvtxTrack* track_out,
698  TVector3& seed_pos,
699  TVector3& seed_mom, TMatrixDSym& seed_cov, const bool do_smearing)
700 {
701  assert(track_out);
702 
703  seed_cov.ResizeTo(6, 6);
704 
705  seed_pos.SetXYZ(0, 0, 0);
706  // reset the seed resolution to the approximate position resolution of the last detector
707  seed_cov[0][0] = .1 * .1;
708  seed_cov[1][1] = .1 * .1;
709  seed_cov[2][2] = 30 * 30;
710  // for (int i = 0; i < 3; i++)
711  // {
712  // seed_cov[i][i] = _phi_resolution * _phi_resolution;
713  // }
714 
715  seed_mom.SetXYZ(0, 0, 10);
716  for (int i = 3; i < 6; i++)
717  {
718  seed_cov[i][i] = 10;
719  }
720 
721  if (particle)
722  {
723  TVector3 True_mom(particle->get_px(), particle->get_py(),
724  particle->get_pz());
725 
726  seed_mom.SetXYZ(particle->get_px(), particle->get_py(),
727  particle->get_pz());
728  if (do_smearing)
729  {
730  const double momSmear = 3. / 180. * M_PI; // rad
731  const double momMagSmear = 0.1; // relative
732 
733  seed_mom.SetMag(
734  True_mom.Mag() + gsl_ran_gaussian(m_RandomGenerator,
735  momMagSmear * True_mom.Mag()));
736  seed_mom.SetTheta(True_mom.Theta() + gsl_ran_gaussian(m_RandomGenerator, momSmear));
737  seed_mom.SetPhi(True_mom.Phi() + gsl_ran_gaussian(m_RandomGenerator, momSmear));
738  }
739  }
740 
741  if (Verbosity())
742  {
743  std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - DEBUG: "
744  << "searching for hits from " << m_PHG4HitContainer.size() << " PHG4Hit nodes" << endl;
745  }
746 
747  // order measurement with g4hit time via stl multimap
748  multimap<double, PHGenFit::Measurement*> ordered_measurements;
749 
750  for (unsigned int ilayer = 0; ilayer < m_PHG4HitContainer.size(); ilayer++)
751  {
752  if (!m_PHG4HitContainer[ilayer])
753  {
754  LogError("No m_PHG4HitContainer[i] found!");
755  continue;
756  }
757 
758  int dettype = m_phg4_detector_type[ilayer];
759  float detradres = m_phg4_detector_radres[ilayer];
760  float detphires = m_phg4_detector_phires[ilayer];
761  float detlonres = m_phg4_detector_lonres[ilayer];
762  float dethiteff = m_phg4_detector_hitfindeff[ilayer];
763  float detnoise = m_phg4_detector_noise[ilayer];
764  if (Verbosity())
765  {
766  std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - DEBUG: "
767  << "ilayer: "
768  << ilayer << ", " << m_PHG4HitsNames[ilayer]
769  << " with nsublayers: " << m_PHG4HitContainer[ilayer]->num_layers()
770  << ", detradres = " << detradres
771  << ", detphires = " << detphires
772  << ", detlonres = " << detlonres
773  << ", dethiteff = " << dethiteff
774  << ", detnoise = " << detnoise
775  << " \n";
776  }
777  for (PHG4HitContainer::LayerIter layerit =
778  m_PHG4HitContainer[ilayer]->getLayers().first;
779  layerit != m_PHG4HitContainer[ilayer]->getLayers().second; layerit++)
780  {
782  m_PHG4HitContainer[ilayer]->getHits(*layerit).first;
783  itr != m_PHG4HitContainer[ilayer]->getHits(*layerit).second; ++itr)
784  {
785  PHG4Hit* hit = itr->second;
786  if (!hit)
787  {
788  LogDebug("No PHG4Hit Found!");
789  continue;
790  }
791 
792  if (hit->get_trkid() == particle->get_track_id() || gsl_ran_binomial(m_RandomGenerator, detnoise, 1) > 0)
793  {
794  if (gsl_ran_binomial(m_RandomGenerator, dethiteff, 1) > 0)
795  {
796  PHGenFit::Measurement* meas = nullptr;
797  if (dettype == Vertical_Plane)
798  {
799  if (Verbosity())
800  {
801  std::cout << "PHG4TrackFastSim::PseudoPatternRecognition -adding vertical plane hit ilayer: "
802  << ilayer << "; detphires: " << detphires << "; detradres: " << detradres << " \n";
803  hit->identify();
804  }
806  detphires, detradres);
807  }
808  else if (dettype == Cylinder)
809  {
810  if (Verbosity())
811  {
812  std::cout << "PHG4TrackFastSim::PseudoPatternRecognition -adding cylinder hit ilayer: "
813  << ilayer << "; detphires: " << detphires << "; detlonres : " << detlonres << " \n";
814  hit->identify();
815  }
816  meas = PHG4HitToMeasurementCylinder(hit,
817  detphires, detlonres);
818  }
819  else
820  {
821  LogError("Type not implemented!");
823  }
824  // meas_out.push_back(meas);
825  ordered_measurements.insert(make_pair(hit->get_avg_t(), meas));
826  track_out->add_g4hit_id(m_PHG4HitContainer[ilayer]->GetID(), hit->get_hit_id());
827  //meas->getMeasurement()->Print(); //DEBUG
828  }
829  }
830  }
831  } /*Loop layers within one detector layer*/
832  } /*Loop detector layers*/
833 
834  for (auto& pair : ordered_measurements)
835  {
836  meas_out.push_back(pair.second);
837 
838  if (Verbosity())
839  {
840  std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - measruement at t = " << pair.first << " ns: ";
841  pair.second->getMeasurement()->Print();
842  }
843  }
844 
845  if (Verbosity())
846  {
847  std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - meas_out.size = " << meas_out.size() << " for "
848  << "particle: "
849  << endl;
850  particle->identify();
851 
852  std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - seed_pos = "
853  << seed_pos.x() << ", " << seed_pos.y() << ". " << seed_pos.z() << endl;
854  std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - seed_pos = "
855  << seed_mom.x() << ", " << seed_mom.y() << ". " << seed_mom.z() << endl;
856  std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - seed_cov = "
857  << sqrt(seed_cov[0][0]) << ", " << sqrt(seed_cov[1][1]) << ". " << sqrt(seed_cov[2][2])
858  << ","
859  << sqrt(seed_cov[3][3]) << ", " << sqrt(seed_cov[4][4]) << ". " << sqrt(seed_cov[5][5]) << endl;
860  }
861 
863 }
864 
866  const PHGenFit::Track* phgf_track,
867  const unsigned int truth_track_id,
868  const unsigned int nmeas,
869  const TVector3& vtx)
870 {
871  assert(out_track);
872 
873  double chi2 = phgf_track->get_chi2();
874  double ndf = phgf_track->get_ndf();
875 
876  double pathlenth_from_first_meas = -999999;
877  unique_ptr<genfit::MeasuredStateOnPlane> gf_state(new genfit::MeasuredStateOnPlane());
878 
879  // if (_detector_type == Vertical_Plane)
880  // {
881  // pathlenth_orig_from_first_meas = phgf_track->extrapolateToPlane(*gf_state, vtx,
882  // TVector3(0., 0., 1.), 0);
883  // }
884  // else if (_detector_type == Cylinder)
885  // pathlenth_orig_from_first_meas = phgf_track->extrapolateToLine(*gf_state, vtx,
886  // TVector3(0., 0., 1.));
887  // else
888  // {
889  // LogError("Detector Type NOT implemented!");
890  // return nullptr;
891  // }
892 
893  // always extrapolate to a z-line through the vertex
894  double pathlenth_orig_from_first_meas = phgf_track->extrapolateToLine(*gf_state, vtx,
895  TVector3(0., 0., 1.));
896 
897  if (Verbosity() > 1)
898  {
899  cout << __PRETTY_FUNCTION__ << __LINE__ << " pathlenth_orig_from_first_meas = " << pathlenth_orig_from_first_meas << endl;
900  }
901  if (pathlenth_orig_from_first_meas < -999990)
902  {
903  LogError("Extraction faild!");
904  return false;
905  }
906 
907  TVector3 mom = gf_state->getMom();
908  TVector3 pos = gf_state->getPos();
909  TMatrixDSym cov = gf_state->get6DCov();
910 
911  out_track->set_truth_track_id(truth_track_id);
918  double dca2d = gf_state->getState()[3];
919  out_track->set_dca2d(dca2d);
920  out_track->set_dca2d_error(gf_state->getCov()[3][3]);
921  double dca3d = sqrt(dca2d * dca2d + gf_state->getState()[4] * gf_state->getState()[4]);
922  out_track->set_dca(dca3d);
923 
924  out_track->set_chisq(chi2);
925  out_track->set_ndf(ndf);
926  out_track->set_charge(phgf_track->get_charge());
927 
928  out_track->set_num_measurements(nmeas);
929 
930  out_track->set_px(mom.Px());
931  out_track->set_py(mom.Py());
932  out_track->set_pz(mom.Pz());
933 
934  out_track->set_x(pos.X());
935  out_track->set_y(pos.Y());
936  out_track->set_z(pos.Z());
937  for (int i = 0; i < 6; i++)
938  {
939  for (int j = i; j < 6; j++)
940  {
941  out_track->set_error(i, j, cov[i][j]);
942  }
943  }
944  // the default name is UNKNOWN - let's set this to ORIGIN since it is at pathlength=0
945  out_track->begin_states()->second->set_name("ORIGIN");
946 
947  // make the projections for all detector types
948  for (map<string, pair<int, double>>::iterator iter = m_ProjectionsMap.begin(); iter != m_ProjectionsMap.end(); ++iter)
949  {
950  switch (iter->second.first)
951  {
952  case DETECTOR_TYPE::Cylinder:
953  pathlenth_from_first_meas = phgf_track->extrapolateToCylinder(*gf_state, iter->second.second, TVector3(0., 0., 0.), TVector3(0., 0., 1.), nmeas - 1);
954  break;
955  case DETECTOR_TYPE::Vertical_Plane:
956  pathlenth_from_first_meas = phgf_track->extrapolateToPlane(*gf_state, TVector3(0., 0., iter->second.second), TVector3(0, 0., 1.), nmeas - 1);
957  break;
958  default:
959  cout << "how in the world did you get here??????" << endl;
960  gSystem->Exit(1);
961  }
962  if (pathlenth_from_first_meas < -999990)
963  {
964  continue;
965  }
966  SvtxTrackState* state = new SvtxTrackState_v1(pathlenth_from_first_meas - pathlenth_orig_from_first_meas);
967  state->set_x(gf_state->getPos().x());
968  state->set_y(gf_state->getPos().y());
969  state->set_z(gf_state->getPos().z());
970 
971  state->set_px(gf_state->getMom().x());
972  state->set_py(gf_state->getMom().y());
973  state->set_pz(gf_state->getMom().z());
974 
975  state->set_name(iter->first);
976  for (int i = 0; i < 6; i++)
977  {
978  for (int j = i; j < 6; j++)
979  {
980  state->set_error(i, j, gf_state->get6DCov()[i][j]);
981  }
982  }
983  out_track->insert_state(state);
984  // the state is cloned on insert_state, so delete this copy here!
985  delete state;
986  }
987 
988  return true;
989 }
990 
992  const PHG4Hit* g4hit, const double phi_resolution,
993  const double r_resolution)
994 {
995  TVector3 pos(g4hit->get_avg_x(), g4hit->get_avg_y(), g4hit->get_avg_z());
996 
997  TVector3 v(pos.X(), pos.Y(), 0);
998  v = 1 / v.Mag() * v;
999 
1000  TVector3 u = v.Cross(TVector3(0, 0, 1));
1001  u = 1 / u.Mag() * u;
1002 
1003  double u_smear = 0.;
1004  double v_smear = 0.;
1005  if (m_SmearingFlag)
1006  {
1007  u_smear = gsl_ran_gaussian(m_RandomGenerator, phi_resolution);
1008  v_smear = gsl_ran_gaussian(m_RandomGenerator, r_resolution);
1009  }
1010  pos.SetX(g4hit->get_avg_x() + u_smear * u.X() + v_smear * v.X());
1011  pos.SetY(g4hit->get_avg_y() + u_smear * u.Y() + v_smear * v.Y());
1012 
1013  PHGenFit::PlanarMeasurement* meas = new PHGenFit::PlanarMeasurement(pos, u, v, phi_resolution,
1014  r_resolution);
1015 
1016  // std::cout<<"------------\n";
1017  // pos.Print();
1018  // std::cout<<"at "<<istation<<" station, "<<ioctant << " octant \n";
1019  // u.Print();
1020  // v.Print();
1021 
1022  //dynamic_cast<PHGenFit::PlanarMeasurement*> (meas)->getMeasurement()->Print();
1023 
1024  return meas;
1025 }
1026 
1028  const PHG4Hit* g4hit, const double phi_resolution,
1029  const double z_resolution)
1030 {
1031  TVector3 pos(g4hit->get_avg_x(), g4hit->get_avg_y(), g4hit->get_avg_z());
1032 
1033  TVector3 v(0, 0, 1);
1034 
1035  TVector3 u = v.Cross(TVector3(pos.X(), pos.Y(), 0));
1036  u = 1 / u.Mag() * u;
1037 
1038  double u_smear = 0.;
1039  double v_smear = 0.;
1040  if (m_SmearingFlag)
1041  {
1042  u_smear = gsl_ran_gaussian(m_RandomGenerator, phi_resolution);
1043  v_smear = gsl_ran_gaussian(m_RandomGenerator, z_resolution);
1044  }
1045  pos.SetX(g4hit->get_avg_x() + u_smear * u.X());
1046  pos.SetY(g4hit->get_avg_y() + u_smear * u.Y());
1047  pos.SetZ(g4hit->get_avg_z() + v_smear);
1048 
1049  PHGenFit::PlanarMeasurement* meas = new PHGenFit::PlanarMeasurement(pos, u, v, phi_resolution,
1050  z_resolution);
1051 
1052  // std::cout<<"------------\n";
1053  // pos.Print();
1054  // std::cout<<"at "<<istation<<" station, "<<ioctant << " octant \n";
1055  // u.Print();
1056  // v.Print();
1057 
1058  //dynamic_cast<PHGenFit::PlanarMeasurement*> (meas)->getMeasurement()->Print();
1059 
1060  return meas;
1061 }
1062 
1063 PHGenFit::Measurement* PHG4TrackFastSim::VertexMeasurement(const TVector3& vtx, double dxy, double dz)
1064 {
1065  TMatrixDSym cov(3);
1066  cov.Zero();
1067  cov(0, 0) = dxy * dxy;
1068  cov(1, 1) = dxy * dxy;
1069  cov(2, 2) = dz * dz;
1070 
1071  TVector3 pos = vtx;
1072  pos.SetX(vtx.X());
1073  pos.SetY(vtx.Y());
1074  pos.SetZ(vtx.Z());
1075 
1077 
1078  return meas;
1079 }
1080 
1082 {
1084  {
1086  }
1087  return;
1088 }
1089 
1090 void PHG4TrackFastSim::add_state_name(const std::string& stateName)
1091 {
1093  {
1094  m_ProjectionsMap.insert(make_pair(stateName, make_pair(DETECTOR_TYPE::Vertical_Plane, NAN)));
1095  }
1097  {
1098  m_ProjectionsMap.insert(make_pair(stateName, make_pair(DETECTOR_TYPE::Cylinder, NAN)));
1099  }
1100  else
1101  {
1102  cout << PHWHERE << " Invalid stateName " << stateName << endl;
1103  cout << endl
1104  << "These are implemented for cylinders" << endl;
1105  for (auto iter : reserved_cylinder_projection_names)
1106  {
1107  cout << iter << endl;
1108  }
1109  cout << endl
1110  << "These are implemented are for zplanes" << endl;
1111  for (auto iter : reserved_zplane_projection_names)
1112  {
1113  cout << iter << endl;
1114  }
1115  gSystem->Exit(1);
1116  }
1117  return;
1118 }
1119 
1120 void PHG4TrackFastSim::add_cylinder_state(const std::string& stateName, const double radius)
1121 {
1124  {
1125  cout << PHWHERE << ": " << stateName << " is a reserved name, used a different name for your cylinder projection" << endl;
1126  gSystem->Exit(1);
1127  }
1128  if (m_ProjectionsMap.find(stateName) != m_ProjectionsMap.end())
1129  {
1130  cout << PHWHERE << ": " << stateName << " is already a projection, please rename" << endl;
1131  gSystem->Exit(1);
1132  }
1133  m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Cylinder, radius)));
1134  return;
1135 }
1136 
1137 void PHG4TrackFastSim::add_zplane_state(const std::string& stateName, const double zplane)
1138 {
1141  {
1142  cout << PHWHERE << ": " << stateName << " is a reserved name, used different name for your zplane projection" << endl;
1143  gSystem->Exit(1);
1144  }
1145  if (m_ProjectionsMap.find(stateName) != m_ProjectionsMap.end())
1146  {
1147  cout << PHWHERE << ": " << stateName << " is already a projection, please rename" << endl;
1148  gSystem->Exit(1);
1149  }
1150  m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Vertical_Plane, zplane)));
1151  return;
1152 }