27 : FW::
BareAlgorithm(
"TruthVerticesToTracksAlgorithm", level), m_cfg(cfg) {
29 throw std::invalid_argument(
"Missing input collection");
31 throw std::invalid_argument(
"Missing output collection");
33 throw std::invalid_argument(
"Missing random number service");
39 const auto& vertexCollection =
40 context.
eventStore.
get<std::vector<FW::SimVertex>>(m_cfg.input);
42 std::shared_ptr<Acts::PerigeeSurface> perigeeSurface =
43 Acts::Surface::makeShared<Acts::PerigeeSurface>(m_cfg.refPosition);
64 std::vector<VertexAndTracks> vertexAndTracksCollection;
67 for (
auto& vtx : vertexCollection) {
71 vertexAndTracks.
vertex = vtx;
74 std::vector<Acts::BoundParameters> trackCollection;
77 for (
auto const&
particle : vtx.outgoing) {
86 auto result = propagator.
propagate(
start, *perigeeSurface, pOptions);
92 const auto& perigeeParameters = (*result).endParameters->parameters();
94 auto newTrackParams = perigeeParameters;
96 if (m_cfg.doSmearing) {
100 m_cfg.ipResA * std::exp(-m_cfg.ipResB * particlePt) + m_cfg.ipResC;
104 std::normal_distribution<double> gaussDist_IP(0., ipRes);
106 std::normal_distribution<double> gaussDist_angular(0., m_cfg.angRes);
108 std::normal_distribution<double> gaussDist_qp(
109 0., m_cfg.qpRelRes * perigeeParameters[4]);
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);
118 smrdParamVec << rn_d0, rn_z0, rn_ph, rn_th, rn_qp, 0.;
121 newTrackParams += smrdParamVec;
124 correctPhiThetaPeriodicity(newTrackParams[2], newTrackParams[3]);
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.;
133 context.
geoContext, covMat, newTrackParams, perigeeSurface));
136 context.
geoContext, std::nullopt, newTrackParams, perigeeSurface));
141 vertexAndTracks.
tracks = trackCollection;
143 vertexAndTracksCollection.push_back(vertexAndTracks);
148 context.
eventStore.
add(m_cfg.output, std::move(vertexAndTracksCollection));
154 double& phiIn,
double& thetaIn)
const {
155 double tmpPhi = std::fmod(phiIn, 2 *
M_PI);
159 if (tmpPhi < -M_PI && tmpPhi > -2 *
M_PI) {
163 double tmpTht = std::fmod(thetaIn, 2 *
M_PI);
164 if (tmpTht < -
M_PI) {
166 }
else if (tmpTht < 0) {
169 tmpPhi = tmpPhi >
M_PI ? tmpPhi - 2 *
M_PI : tmpPhi;
172 tmpTht = 2 *
M_PI - tmpTht;
174 tmpPhi = tmpPhi >
M_PI ? (tmpPhi - 2 *
M_PI) : tmpPhi;