ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHActsTrackProjection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHActsTrackProjection.cc
2 
5 #include <phool/getClass.h>
6 #include <phool/phool.h>
7 #include <phool/PHDataNode.h>
8 #include <phool/PHNode.h>
9 #include <phool/PHNodeIterator.h>
10 #include <phool/PHObject.h>
11 #include <phool/PHTimer.h>
12 
17 
18 #include <calobase/RawTowerGeomContainer.h>
19 #include <calobase/RawTowerContainer.h>
20 #include <calobase/RawTower.h>
21 #include <calobase/RawClusterContainer.h>
22 #include <calobase/RawCluster.h>
23 #include <calobase/RawClusterUtility.h>
24 #include <phgeom/PHGeomUtility.h>
25 
26 #include <Acts/Geometry/GeometryIdentifier.hpp>
32 
33 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
34 
36 #include <math.h>
37 
39  : SubsysReco(name)
40  , m_trajectories(nullptr)
41 {
42  m_caloNames.push_back("CEMC");
43  m_caloNames.push_back("HCALIN");
44  m_caloNames.push_back("HCALOUT");
45 
46  m_caloTypes.push_back(SvtxTrack::CEMC);
47  m_caloTypes.push_back(SvtxTrack::HCALIN);
49 }
50 
52 {
53  if(Verbosity() > 1)
54  std::cout << "PHActsTrackProjection begin Init" << std::endl;
55 
56  int ret = makeCaloSurfacePtrs(topNode);
57 
60 
61  if(Verbosity() > 1)
62  std::cout << "PHActsTrackProjection finished Init" << std::endl;
63 
64  return ret;
65 }
66 
68 {
69  if(Verbosity() > 1)
70  std::cout << "PHActsTrackProjection : Starting process_event event "
71  << m_event << std::endl;
72 
73  for(int layer = 0; layer < m_nCaloLayers; layer++)
74  {
75  if(Verbosity() > 1)
76  std::cout << "Processing calo layer "
77  << m_caloNames.at(layer) << std::endl;
78 
79  int ret = projectTracks(topNode, layer);
82  }
83 
84  if(Verbosity() > 1)
85  std::cout << "PHActsTrackProjection : Finished process_event event "
86  << m_event << std::endl;
87 
88  m_event++;
89 
91 }
92 
94 {
96 }
97 
99 {
101 }
102 
103 
105  const int caloLayer)
106 {
107  if (setCaloContainerNodes(topNode, caloLayer)
110 
111  for(const auto& [trackKey, traj] : *m_trajectories)
112  {
113  const auto track = m_trackMap->get(trackKey);;
114  const auto& [trackTips, mj] = traj.trajectory();
115  if(trackTips.size() > 1 and Verbosity() > 0)
116  {
117  std::cout << PHWHERE
118  << "More than 1 track tip per track. Should never happen..."
119  << std::endl;
120  }
121  for(const auto& trackTip : trackTips)
122  {
123  const auto params = traj.trackParameters(trackTip);
124  auto cylSurf =
125  m_caloSurfaces.find(m_caloNames.at(caloLayer))->second;
126 
127  auto result = propagateTrack(params, cylSurf);
128 
129  if(result.ok())
130  {
131  auto trackStateParams = std::move(**result);
132  updateSvtxTrack(trackStateParams,
133  track, caloLayer);
134  }
135  }
136  }
137 
139 }
140 
143 {
144 
145  Acts::Vector3D momentum(track->get_px(),
146  track->get_py(),
147  track->get_pz());
148 
149  auto actsVertex = getVertex(track);
150  auto perigee =
151  Acts::Surface::makeShared<Acts::PerigeeSurface>(actsVertex);
152  auto actsFourPos =
153  Acts::Vector4D(track->get_x() * Acts::UnitConstants::cm,
154  track->get_y() * Acts::UnitConstants::cm,
155  track->get_z() * Acts::UnitConstants::cm,
157 
159  for(int i = 0; i < 6; i++)
160  for(int j = 0; j < 6; j++)
161  cov(i,j) = track->get_acts_covariance(i,j);
162 
164  actsFourPos, momentum,
165  track->get_p(), track->get_charge(),
166  cov);
167  return param;
168 }
170 {
171  auto vertexId = track->get_vertex_id();
172  const SvtxVertex* svtxVertex = m_vertexMap->get(vertexId);
173 
174  Acts::Vector3D vertex(svtxVertex->get_x() * Acts::UnitConstants::cm,
175  svtxVertex->get_y() * Acts::UnitConstants::cm,
176  svtxVertex->get_z() * Acts::UnitConstants::cm);
177  return vertex;
178 }
179 
180 
182  const Acts::BoundTrackParameters& params,
183  SvtxTrack* svtxTrack,
184  const int caloLayer)
185 {
186  auto projectionPos = params.position(m_tGeometry->geoContext);
187 
188  auto projectionPhi = atan2(projectionPos(1), projectionPos(0));
189  auto projectionEta = asinh(projectionPos(2) /
190  sqrt(projectionPos(0) * projectionPos(0)
191  + projectionPos(1) * projectionPos(1)));
192 
193 
194  if(Verbosity() > 2)
195  {
196  double radius = sqrt(projectionPos(0) * projectionPos(0)
197  + projectionPos(1) * projectionPos(1));
198  std::cout << "Track projection phi/eta " << projectionPhi
199  << " and " << projectionEta << std::endl
200  << " projection position " << projectionPos.transpose()
201  << " and radius is " << radius << std::endl;
202  }
203 
204  if(fabs(projectionEta) >= 1.1)
205  return;
206 
207  auto phiBin = m_towerGeomContainer->get_phibin(projectionPhi);
208  auto etaBin = m_towerGeomContainer->get_etabin(projectionEta);
209 
210  double energy3x3 = 0.0;
211  double energy5x5 = 0.0;
212 
213  getSquareTowerEnergies(phiBin, etaBin, energy3x3, energy5x5);
214 
215  if(Verbosity() > 2)
216  std::cout << "3x3/5x5 energy sums " << energy3x3
217  << " and " << energy5x5 << std::endl;
218 
219  svtxTrack->set_cal_energy_3x3(m_caloTypes.at(caloLayer),
220  energy3x3);
221  svtxTrack->set_cal_energy_5x5(m_caloTypes.at(caloLayer),
222  energy5x5);
223 
224  double minIndex = NAN;
225  double minDphi = NAN;
226  double minDeta = NAN;
227  double minE = NAN;
228 
229  getClusterProperties(projectionPhi, projectionEta,
230  minIndex, minDphi, minDeta, minE);
231 
232  if(!std::isnan(minIndex))
233  {
234  svtxTrack->set_cal_dphi(m_caloTypes.at(caloLayer),
235  minDphi);
236  svtxTrack->set_cal_deta(m_caloTypes.at(caloLayer),
237  minDeta);
238  svtxTrack->set_cal_cluster_id(m_caloTypes.at(caloLayer),
239  minIndex);
240  svtxTrack->set_cal_cluster_e(m_caloTypes.at(caloLayer),
241  minE);
242  if(Verbosity() > 2)
243  std::cout << "Calo cluster has dphi " << minDphi
244  << " and deta " << minDeta << " and clusID "
245  << minIndex << " and energy " << minE
246  << std::endl;
247 
248  }
249 
250  return;
251 }
252 
254  double eta,
255  double& minIndex,
256  double& minDphi,
257  double& minDeta,
258  double& minE)
259 {
260  double minR = DBL_MAX;
261  auto clusterMap = m_clusterContainer->getClustersMap();
262  for(const auto& [key, cluster] : clusterMap)
263  {
264  const auto clusterEta =
266  CLHEP::Hep3Vector(0,0,0));
267  const auto dphi = deltaPhi(phi - cluster->get_phi());
268  const auto deta = eta - clusterEta;
269  const auto r = sqrt(pow(dphi,2) + pow(deta,2));
270 
271  if(r < minR)
272  {
273  minIndex = key;
274  minR = r;
275  minDphi = dphi;
276  minDeta = deta;
277  minE = cluster->get_energy();
278  }
279  }
280 
281  return;
282 }
283 
285  int etaBin,
286  double& energy3x3,
287  double& energy5x5)
288 {
289  for(int iphi = phiBin - 2; iphi <= phiBin + 2; ++iphi)
290  {
291  for(int ieta = etaBin - 2; ieta <= etaBin + 2; ++ieta)
292  {
294  int wrapPhi = iphi;
295  if(wrapPhi < 0)
296  wrapPhi += m_towerGeomContainer->get_phibins();
297  if(wrapPhi >= m_towerGeomContainer->get_phibins())
298  wrapPhi -= m_towerGeomContainer->get_phibins();
299 
301  if (ieta < 0 or ieta >= m_towerGeomContainer->get_etabins())
302  continue;
303 
304  auto tower = m_towerContainer->getTower(ieta, wrapPhi);
305 
306  if(!tower)
307  continue;
308 
309  energy5x5 += tower->get_energy();
310  if(abs(iphi - phiBin) <= 1 and abs(ieta - etaBin) <= 1)
311  energy3x3 += tower->get_energy();
312 
313  }
314  }
315 
316  return;
317 }
318 
320  const Acts::BoundTrackParameters& params,
321  const SurfacePtr& targetSurf)
322 {
323 
324  if(Verbosity() > 1)
325  std::cout << "Propagating final track fit with momentum: "
326  << params.momentum() << " and position "
327  << params.position(m_tGeometry->geoContext)
328  << std::endl
329  << "track fit phi/eta "
330  << atan2(params.momentum()(1),
331  params.momentum()(0))
332  << " and "
333  << atanh(params.momentum()(2)
334  / params.momentum().norm())
335  << std::endl;
336 
337  return std::visit([params, targetSurf, this]
338  (auto && inputField) -> BoundTrackParamPtrResult {
339  using InputMagneticField =
340  typename std::decay_t<decltype(inputField)>::element_type;
341  using MagneticField = Acts::SharedBField<InputMagneticField>;
344 
345  MagneticField field(inputField);
346  Stepper stepper(field);
347  Propagator propagator(stepper);
348 
350  if(Verbosity() > 3)
351  logLevel = Acts::Logging::VERBOSE;
352 
353  auto logger = Acts::getDefaultLogger("PHActsTrackProjection",
354  logLevel);
355 
358  Acts::LoggerWrapper{*logger});
359 
360  auto result = propagator.propagate(params, *targetSurf,
361  options);
362 
363  if(result.ok())
364  return std::move((*result).endParameters);
365 
366  return result.error();
367  },
368  std::move(m_tGeometry->magField));
369 
370 }
371 
373  const int caloLayer)
374 {
375  std::string towerGeoNodeName = "TOWERGEOM_" + m_caloNames.at(caloLayer);
376  std::string towerNodeName = "TOWER_CALIB_" + m_caloNames.at(caloLayer);
377  std::string clusterNodeName = "CLUSTER_" + m_caloNames.at(caloLayer);
378 
379  m_towerGeomContainer = findNode::getClass<RawTowerGeomContainer>
380  (topNode, towerGeoNodeName.c_str());
381 
382  m_towerContainer = findNode::getClass<RawTowerContainer>
383  (topNode, towerNodeName.c_str());
384 
385  m_clusterContainer = findNode::getClass<RawClusterContainer>
386  (topNode, clusterNodeName.c_str());
387 
388  if(m_useCemcPosRecalib and
389  m_caloNames.at(caloLayer).compare("CEMC") == 0)
390  {
391  std::string nodeName = "CLUSTER_POS_COR_" + m_caloNames.at(caloLayer);
392  m_clusterContainer = findNode::getClass<RawClusterContainer>
393  (topNode, nodeName.c_str());
394  }
395 
396  if(!m_towerGeomContainer or !m_towerContainer or !m_clusterContainer)
397  {
398  std::cout << PHWHERE
399  << "Calo geometry and/or cluster container not found on node tree. Bailing."
400  << std::endl;
402  }
403 
405 }
406 
408 {
409  for(int caloLayer = 0; caloLayer < m_nCaloLayers; caloLayer++)
410  {
411  if(setCaloContainerNodes(topNode, caloLayer) != Fun4AllReturnCodes::EVENT_OK)
413 
414  const auto caloRadius = m_towerGeomContainer->get_radius()
418  const auto eta = 2.5;
419  const auto theta = 2. * atan(exp(-eta));
420  const auto halfZ = caloRadius / tan(theta)
422 
424  auto transform = Acts::Transform3D::Identity();
425 
426  std::shared_ptr<Acts::CylinderSurface> surf =
427  Acts::Surface::makeShared<Acts::CylinderSurface>(transform,
428  caloRadius,
429  halfZ);
430 
431  m_caloSurfaces.insert(std::make_pair(m_caloNames.at(caloLayer),
432  surf));
433  }
434 
435  if(Verbosity() > 1)
436  {
437  for(const auto& [name, surfPtr] : m_caloSurfaces)
438  {
439  std::cout << "Cylinder " << name << " has center "
440  << surfPtr.get()->center(m_tGeometry->geoContext)
441  << std::endl;
442  }
443  }
444 
446 
447 }
448 
450 {
451  m_trajectories = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode, "ActsTrajectories");
452  if(!m_trajectories)
453  {
454  std::cout << PHWHERE << "No Acts trajectories on node tree, bailing."
455  << std::endl;
457  }
458 
459  m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
460  if(!m_vertexMap)
461  {
462  std::cout << PHWHERE << "No vertex map on node tree, bailing."
463  << std::endl;
465  }
466 
467  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(
468  topNode, "ActsTrackingGeometry");
469  if(!m_tGeometry)
470  {
471  std::cout << "ActsTrackingGeometry not on node tree. Exiting."
472  << std::endl;
473 
475  }
476 
477  m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
478  if(!m_trackMap)
479  {
480  std::cout << PHWHERE << "No SvtxTrackMap on node tree. Bailing."
481  << std::endl;
483  }
484 
486 }
487 
489 {
490  if (phi > M_PI)
491  return phi - 2. * M_PI;
492  else if (phi <= -M_PI)
493  return phi + 2.* M_PI;
494  else
495  return phi;
496 }