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
1 #include "PHActsTrackProjection.h"
2 
3 #include <fun4all/Fun4AllReturnCodes.h>
4 #include <phool/PHCompositeNode.h>
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 
13 #include <trackbase_historic/SvtxTrackMap.h>
14 #include <trackbase_historic/SvtxTrackState.h>
15 #include <trackbase_historic/SvtxVertexMap.h>
16 #include <trackbase_historic/SvtxVertex.h>
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>
27 #include <Acts/MagneticField/ConstantBField.hpp>
28 #include <Acts/MagneticField/InterpolatedBFieldMap.hpp>
29 #include <Acts/MagneticField/SharedBField.hpp>
30 #include <Acts/Propagator/EigenStepper.hpp>
31 #include <Acts/Surfaces/PerigeeSurface.hpp>
32 
33 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
34 
35 #include <CLHEP/Vector/ThreeVector.h>
36 #include <math.h>
37 
38 PHActsTrackProjection::PHActsTrackProjection(const std::string& name)
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);
48  m_caloTypes.push_back(SvtxTrack::HCALOUT);
49 }
50 
51 int PHActsTrackProjection::InitRun(PHCompositeNode *topNode)
52 {
53  if(Verbosity() > 1)
54  std::cout << "PHActsTrackProjection begin Init" << std::endl;
55 
56  int ret = makeCaloSurfacePtrs(topNode);
57 
58  if(getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
59  ret = Fun4AllReturnCodes::ABORTEVENT;
60 
61  if(Verbosity() > 1)
62  std::cout << "PHActsTrackProjection finished Init" << std::endl;
63 
64  return ret;
65 }
66 
67 int PHActsTrackProjection::process_event(PHCompositeNode *topNode)
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);
80  if(ret != Fun4AllReturnCodes::EVENT_OK)
81  return Fun4AllReturnCodes::ABORTEVENT;
82  }
83 
84  if(Verbosity() > 1)
85  std::cout << "PHActsTrackProjection : Finished process_event event "
86  << m_event << std::endl;
87 
88  m_event++;
89 
90  return Fun4AllReturnCodes::EVENT_OK;
91 }
92 
93 int PHActsTrackProjection::Init(PHCompositeNode */*topNode*/)
94 {
95  return Fun4AllReturnCodes::EVENT_OK;
96 }
97 
98 int PHActsTrackProjection::End(PHCompositeNode */*topNode*/)
99 {
100  return Fun4AllReturnCodes::EVENT_OK;
101 }
102 
103 
104 int PHActsTrackProjection::projectTracks(PHCompositeNode *topNode,
105  const int caloLayer)
106 {
107  if (setCaloContainerNodes(topNode, caloLayer)
108  != Fun4AllReturnCodes::EVENT_OK)
109  return Fun4AllReturnCodes::ABORTEVENT;
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 
138  return Fun4AllReturnCodes::EVENT_OK;
139 }
140 
141 Acts::BoundTrackParameters
142 PHActsTrackProjection::makeTrackParams(SvtxTrack* track)
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,
156  10 * Acts::UnitConstants::ns);
157 
158  Acts::BoundSymMatrix cov;
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 
163  Acts::BoundTrackParameters param(perigee, m_tGeometry->geoContext,
164  actsFourPos, momentum,
165  track->get_p(), track->get_charge(),
166  cov);
167  return param;
168 }
169 Acts::Vector3D PHActsTrackProjection::getVertex(SvtxTrack *track)
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 
181 void PHActsTrackProjection::updateSvtxTrack(
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 
253 void PHActsTrackProjection::getClusterProperties(double phi,
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 =
265  RawClusterUtility::GetPseudorapidity(*cluster,
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 
284 void PHActsTrackProjection::getSquareTowerEnergies(int phiBin,
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 
319 BoundTrackParamPtrResult PHActsTrackProjection::propagateTrack(
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>;
342  using Stepper = Acts::EigenStepper<MagneticField>;
343  using Propagator = Acts::Propagator<Stepper>;
344 
345  MagneticField field(inputField);
346  Stepper stepper(field);
347  Propagator propagator(stepper);
348 
349  Acts::Logging::Level logLevel = Acts::Logging::FATAL;
350  if(Verbosity() > 3)
351  logLevel = Acts::Logging::VERBOSE;
352 
353  auto logger = Acts::getDefaultLogger("PHActsTrackProjection",
354  logLevel);
355 
356  Acts::PropagatorOptions<> options(m_tGeometry->geoContext,
357  m_tGeometry->magFieldContext,
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 
372 int PHActsTrackProjection::setCaloContainerNodes(PHCompositeNode *topNode,
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;
401  return Fun4AllReturnCodes::ABORTEVENT;
402  }
403 
404  return Fun4AllReturnCodes::EVENT_OK;
405 }
406 
407 int PHActsTrackProjection::makeCaloSurfacePtrs(PHCompositeNode *topNode)
408 {
409  for(int caloLayer = 0; caloLayer < m_nCaloLayers; caloLayer++)
410  {
411  if(setCaloContainerNodes(topNode, caloLayer) != Fun4AllReturnCodes::EVENT_OK)
412  return Fun4AllReturnCodes::ABORTEVENT;
413 
414  const auto caloRadius = m_towerGeomContainer->get_radius()
415  * Acts::UnitConstants::cm;
418  const auto eta = 2.5;
419  const auto theta = 2. * atan(exp(-eta));
420  const auto halfZ = caloRadius / tan(theta)
421  * Acts::UnitConstants::cm;
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 
445  return Fun4AllReturnCodes::EVENT_OK;
446 
447 }
448 
449 int PHActsTrackProjection::getNodes(PHCompositeNode *topNode)
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;
456  return Fun4AllReturnCodes::ABORTEVENT;
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;
464  return Fun4AllReturnCodes::ABORTEVENT;
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 
474  return Fun4AllReturnCodes::ABORTEVENT;
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;
482  return Fun4AllReturnCodes::ABORTEVENT;
483  }
484 
485  return Fun4AllReturnCodes::EVENT_OK;
486 }
487 
488 double PHActsTrackProjection::deltaPhi(const double& phi)
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 }