ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TruthVerticesToTracks.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TruthVerticesToTracks.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
10 
11 #include <iostream>
12 #include <optional>
13 #include <stdexcept>
14 
23 
27  : FW::BareAlgorithm("TruthVerticesToTracksAlgorithm", level), m_cfg(cfg) {
28  if (m_cfg.input.empty()) {
29  throw std::invalid_argument("Missing input collection");
30  } else if (m_cfg.output.empty()) {
31  throw std::invalid_argument("Missing output collection");
32  } else if (m_cfg.randomNumberSvc == nullptr) {
33  throw std::invalid_argument("Missing random number service");
34  }
35 }
36 
38  const AlgorithmContext& context) const {
39  const auto& vertexCollection =
40  context.eventStore.get<std::vector<FW::SimVertex>>(m_cfg.input);
41 
42  std::shared_ptr<Acts::PerigeeSurface> perigeeSurface =
43  Acts::Surface::makeShared<Acts::PerigeeSurface>(m_cfg.refPosition);
44 
45  // Set up constant B-Field
46  Acts::ConstantBField bField(m_cfg.bField);
47 
48  // Set up stepper
50 
51  // Set up propagator with void navigator
53  stepper);
54 
55  // Set up propagator options
56  Acts::PropagatorOptions<> pOptions(context.geoContext,
57  context.magFieldContext);
58  pOptions.direction = Acts::backward;
59 
60  // Create random number generator and spawn gaussian distribution
61  FW::RandomEngine rng = m_cfg.randomNumberSvc->spawnGenerator(context);
62 
63  // Vector to store VertexAndTracks extracted from event
64  std::vector<VertexAndTracks> vertexAndTracksCollection;
65 
66  // Start looping over all vertices in current event
67  for (auto& vtx : vertexCollection) {
68  // Create VertexAndTracks object
69  VertexAndTracks vertexAndTracks;
70  // Store current vertex
71  vertexAndTracks.vertex = vtx;
72 
73  // Track objects at current vertex
74  std::vector<Acts::BoundParameters> trackCollection;
75 
76  // Iterate over all particle emerging from current vertex
77  for (auto const& particle : vtx.outgoing) {
78  const Acts::Vector3D& ptclMom =
79  particle.absMomentum() * particle.unitDirection();
80 
81  // Define start track params
82  Acts::CurvilinearParameters start(std::nullopt, particle.position(),
83  ptclMom, particle.charge(),
84  particle.time());
85  // Run propagator
86  auto result = propagator.propagate(start, *perigeeSurface, pOptions);
87  if (!result.ok()) {
88  continue;
89  }
90 
91  // get perigee parameters
92  const auto& perigeeParameters = (*result).endParameters->parameters();
93 
94  auto newTrackParams = perigeeParameters;
95 
96  if (m_cfg.doSmearing) {
97  // Calculate pt-dependent IP resolution
98  const double particlePt = Acts::VectorHelpers::perp(ptclMom);
99  const double ipRes =
100  m_cfg.ipResA * std::exp(-m_cfg.ipResB * particlePt) + m_cfg.ipResC;
101 
102  // except for IP resolution, following variances are rough guesses
103  // Gaussian distribution for IP resolution
104  std::normal_distribution<double> gaussDist_IP(0., ipRes);
105  // Gaussian distribution for angular resolution
106  std::normal_distribution<double> gaussDist_angular(0., m_cfg.angRes);
107  // Gaussian distribution for q/p (momentum) resolution
108  std::normal_distribution<double> gaussDist_qp(
109  0., m_cfg.qpRelRes * perigeeParameters[4]);
110 
111  double rn_d0 = gaussDist_IP(rng);
112  double rn_z0 = gaussDist_IP(rng);
113  double rn_ph = gaussDist_angular(rng);
114  double rn_th = gaussDist_angular(rng);
115  double rn_qp = gaussDist_qp(rng);
116 
117  Acts::BoundVector smrdParamVec;
118  smrdParamVec << rn_d0, rn_z0, rn_ph, rn_th, rn_qp, 0.;
119 
120  // Update track parameters
121  newTrackParams += smrdParamVec;
122 
123  // Correct for phi and theta wrap
124  correctPhiThetaPeriodicity(newTrackParams[2], newTrackParams[3]);
125 
126  // Update track covariance
127  Acts::BoundSymMatrix covMat;
128  covMat.setZero();
129  covMat.diagonal() << rn_d0 * rn_d0, rn_z0 * rn_z0, rn_ph * rn_ph,
130  rn_th * rn_th, rn_qp * rn_qp, 1.;
131 
132  trackCollection.push_back(Acts::BoundParameters(
133  context.geoContext, covMat, newTrackParams, perigeeSurface));
134  } else {
135  trackCollection.push_back(Acts::BoundParameters(
136  context.geoContext, std::nullopt, newTrackParams, perigeeSurface));
137  }
138  } // end iteration over all particle at vertex
139 
140  // Store track objects in VertexAndTracks
141  vertexAndTracks.tracks = trackCollection;
142  // Add to collection
143  vertexAndTracksCollection.push_back(vertexAndTracks);
144 
145  } // end iteration over all vertices
146 
147  // VertexAndTracks objects to the EventStore
148  context.eventStore.add(m_cfg.output, std::move(vertexAndTracksCollection));
149 
151 }
152 
154  double& phiIn, double& thetaIn) const {
155  double tmpPhi = std::fmod(phiIn, 2 * M_PI); // temp phi
156  if (tmpPhi > M_PI) {
157  tmpPhi -= 2 * M_PI;
158  }
159  if (tmpPhi < -M_PI && tmpPhi > -2 * M_PI) {
160  tmpPhi += 2 * M_PI;
161  }
162 
163  double tmpTht = std::fmod(thetaIn, 2 * M_PI); // temp theta
164  if (tmpTht < -M_PI) {
165  tmpTht = std::abs(tmpTht + 2 * M_PI);
166  } else if (tmpTht < 0) {
167  tmpTht *= -1;
168  tmpPhi += M_PI;
169  tmpPhi = tmpPhi > M_PI ? tmpPhi - 2 * M_PI : tmpPhi;
170  }
171  if (tmpTht > M_PI) {
172  tmpTht = 2 * M_PI - tmpTht;
173  tmpPhi += M_PI;
174  tmpPhi = tmpPhi > M_PI ? (tmpPhi - 2 * M_PI) : tmpPhi;
175  }
176 
177  phiIn = tmpPhi;
178  thetaIn = tmpTht;
179 }