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