Run the processor.
151 {
152
153
154 auto tg{geometry()};
155
156
157 if (!event.
exists(trackCollection_))
return;
158 auto tracks{
event.getCollection<
ldmx::Track>(trackCollection_)};
159
160
161 if (!event.
exists(measCollection_))
return;
163
165
166
167 Acts::GainMatrixUpdater updater;
168 Acts::GsfExtensions<Acts::VectorMultiTrajectory> gsf_extensions;
169 gsf_extensions.updater.connect<
170 &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(
171 &updater);
172 gsf_extensions.calibrator
174 Acts::VectorMultiTrajectory>>(&calibrator);
175
176
177 struct SurfaceAccessor {
178 const Acts::TrackingGeometry* trackingGeometry;
179
180 const Acts::Surface* operator()(const Acts::SourceLink& sourceLink) const {
181 const auto& indexSourceLink =
183 return trackingGeometry->findSurface(indexSourceLink.geometryId());
184 }
185 };
186
187 SurfaceAccessor m_slSurfaceAccessor{tg.getTG().get()};
188
189 gsf_extensions.surfaceAccessor.connect<&SurfaceAccessor::operator()>(
190 &m_slSurfaceAccessor);
191 gsf_extensions.mixtureReducer.connect<&Acts::reduceMixtureLargestWeights>();
192
193
194
195
196 Acts::PropagatorOptions<Acts::StepperPlainOptions,
197 Acts::NavigatorPlainOptions, ActionList, AbortList>
198 propagator_options(geometry_context(), magnetic_field_context());
199
200 propagator_options.pathLimit = std::numeric_limits<double>::max();
201
202
203 propagator_options.loopProtection =
204 false;
205
206
207 auto& mInteractor =
208 propagator_options.actionList.get<Acts::MaterialInteractor>();
209 mInteractor.multipleScattering = true;
210 mInteractor.energyLoss = true;
211 mInteractor.recordInteractions = false;
212
213
214 auto& sLogger =
215 propagator_options.actionList.get<Acts::detail::SteppingLogger>();
216 sLogger.sterile = true;
217
218 propagator_options.stepping.maxStepSize =
219 propagator_step_size_ * Acts::UnitConstants::mm;
220 propagator_options.maxSteps = propagator_maxSteps_;
221
222
223
224
225
226
227 std::shared_ptr<const Acts::PerigeeSurface> origin_surface =
228 Acts::Surface::makeShared<Acts::PerigeeSurface>(
229 Acts::Vector3(0., 0., 0.));
230
231 std::shared_ptr<const Acts::PerigeeSurface> tagger_layer_surface =
232 Acts::Surface::makeShared<Acts::PerigeeSurface>(
233 Acts::Vector3(-700., 0., 0.));
234
235 std::shared_ptr<const Acts::PerigeeSurface> reference_surface =
236 origin_surface;
237 if (taggerTracking_) {
238 reference_surface = tagger_layer_surface;
239 }
240
241
242
243
244
245
246
247
248
249 Acts::GsfOptions<Acts::VectorMultiTrajectory> gsfOptions{
250 geometry_context(), magnetic_field_context(), calibration_context()};
251
252 gsfOptions.extensions = gsf_extensions;
253 gsfOptions.propagatorPlainOptions = propagator_options;
254 gsfOptions.referenceSurface = &(*reference_surface);
255 gsfOptions.maxComponents = maxComponents_;
256 gsfOptions.weightCutoff = weightCutoff_;
257 gsfOptions.abortOnError = abortOnError_;
258 gsfOptions.disableAllMaterialHandling = disableAllMaterialHandling_;
259
260
261 std::vector<ldmx::Track> out_tracks;
262
263
264 Acts::VectorTrackContainer vtc;
265 Acts::VectorMultiTrajectory mtj;
266 Acts::TrackContainer tc{vtc, mtj};
267
268
269 unsigned int itrk = 0;
270
271 for (auto& track : tracks) {
272
273 std::vector<ldmx::Measurement> measOnTrack;
274
275
276 std::vector<Acts::SourceLink> fit_trackSourceLinks;
277
278 for (auto imeas : track.getMeasurementsIdxs()) {
279 auto meas = measurements.at(imeas);
280 measOnTrack.push_back(meas);
281
282
283
284 const Acts::Surface* hit_surface = tg.getSurface(meas.getLayerID());
285
286
288 fit_trackSourceLinks.push_back(Acts::SourceLink(idx_sl));
289 }
290
291
292 std::reverse(measOnTrack.begin(), measOnTrack.end());
293 std::reverse(fit_trackSourceLinks.begin(), fit_trackSourceLinks.end());
294
295 for (auto m : measOnTrack) {
296 ldmx_log(debug) << "Measurement:\n" << m << "\n";
297 }
298
299 ldmx_log(debug) << "GSF Refitting";
300
301
302
303 std::shared_ptr<Acts::PerigeeSurface> perigee =
304 Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(
305 track.getPerigeeX(), track.getPerigeeY(), track.getPerigeeZ()));
306
307 Acts::BoundTrackParameters trk_btp =
308 tracking::sim::utils::boundTrackParameters(track, perigee);
309
310 std::shared_ptr<Acts::Surface> beamOrigin_surface =
311 tracking::sim::utils::unboundSurface(-700);
312
313 const std::shared_ptr<Acts::Surface> target_surface =
314 tracking::sim::utils::unboundSurface(0.);
315
316 const std::shared_ptr<Acts::Surface> ecal_surface =
317 tracking::sim::utils::unboundSurface(240.5);
318
319 Acts::BoundTrackParameters trk_btp_bO =
320 tracking::sim::utils::boundTrackParameters(track, perigee);
321
322 if (taggerTracking_) {
323 if (!track.getTrackState(ldmx::TrackStateType::AtBeamOrigin)
324 .has_value()) {
325 ldmx_log(warn) << "Failed retreiving AtBeamOrigin TrackState for "
326 "track. Skipping..";
327 continue;
328 }
329
330 auto ts = track.getTrackState(ldmx::TrackStateType::AtBeamOrigin).value();
331 trk_btp_bO = tracking::sim::utils::btp(
332 ts, beamOrigin_surface,
333 11);
334 } else {
335 if (!track.getTrackState(ldmx::TrackStateType::AtTarget).has_value()) {
336 ldmx_log(warn)
337 << "Failed retreiving AtTarget TrackState for track. Skipping..";
338 continue;
339 }
340 auto ts = track.getTrackState(ldmx::TrackStateType::AtTarget).value();
341 trk_btp_bO = tracking::sim::utils::btp(ts, target_surface, 11);
342 }
343 const Acts::BoundVector& trkpars = trk_btp.parameters();
344 ldmx_log(debug) << "CKF Track parameters" << std::endl
345 << trkpars[0] << " " << trkpars[1] << " " << trkpars[2]
346 << " " << trkpars[3] << " " << trkpars[4] << " "
347 << trkpars[5] << std::endl
348 << "Perigee Surface" << std::endl
349 << track.getPerigeeX() << " " << track.getPerigeeY() << " "
350 << track.getPerigeeZ();
351
352 Acts::Vector3 trk_pos = trk_btp.position(geometry_context());
353
354 ldmx_log(debug) << trk_pos(0) << " " << trk_pos(1) << " " << trk_pos(2)
355 << std::endl;
356
357 const Acts::BoundVector& trkpars_bO = trk_btp_bO.parameters();
358 ldmx_log(debug) << "CKF BeamOrigin track parameters" << std::endl
359 << trkpars_bO[0] << " " << trkpars_bO[1] << " "
360 << trkpars_bO[2] << " " << trkpars_bO[3] << " "
361 << trkpars_bO[4] << " " << trkpars_bO[5] << " ";
362
363 Acts::Vector3 trk_pos_bO = trk_btp_bO.position(geometry_context());
364 ldmx_log(debug) << trk_pos_bO(0) << " " << trk_pos_bO(1) << " "
365 << trk_pos_bO(2) << std::endl;
366
367 auto gsf_refit_result =
368 gsf_->fit(fit_trackSourceLinks.begin(), fit_trackSourceLinks.end(),
369 trk_btp_bO, gsfOptions, tc);
370
371 if (!gsf_refit_result.ok()) {
372 ldmx_log(warn) << "GSF re-fit failed" << std::endl;
373 continue;
374 }
375
376 if (tc.size() < 1) continue;
377
378 auto gsftrk = tc.getTrack(itrk);
379 calculateTrackQuantities(gsftrk);
380
381 const Acts::BoundVector& perigee_pars = gsftrk.parameters();
382 const Acts::BoundMatrix& trk_cov = gsftrk.covariance();
383 const Acts::Surface& perigee_surface = gsftrk.referenceSurface();
384
385 ldmx_log(debug)
386 << "Found track:" << std::endl
387 << "Track states " << gsftrk.nTrackStates() << std::endl
388 << perigee_pars[Acts::eBoundLoc0] << " "
389 << perigee_pars[Acts::eBoundLoc1] << " "
390 << perigee_pars[Acts::eBoundPhi] << " "
391 << perigee_pars[Acts::eBoundTheta] << " "
392 << perigee_pars[Acts::eBoundQOverP] << std::endl
393 << "Reference Surface" << std::endl
394 << " " << perigee_surface.transform(geometry_context()).translation()(0)
395 << " " << perigee_surface.transform(geometry_context()).translation()(1)
396 << " " << perigee_surface.transform(geometry_context()).translation()(2)
397 << std::endl;
398
400
401 bool success = false;
402 if (taggerTracking_) {
403 ldmx_log(debug) << "Target extrapolation";
405
406 success = trk_extrap_->TrackStateAtSurface(
407 gsftrk, target_surface, tsAtTarget, ldmx::TrackStateType::AtTarget);
408
409 if (success) trk.addTrackState(tsAtTarget);
410 } else {
411 ldmx_log(debug) << "Ecal Extrapolation";
413 success = trk_extrap_->TrackStateAtSurface(gsftrk, ecal_surface, tsAtEcal,
414 ldmx::TrackStateType::AtECAL);
415
416 if (success) trk.addTrackState(tsAtEcal);
417 }
418
419 trk.setPerigeeLocation(
420 perigee_surface.transform(geometry_context()).translation()(0),
421 perigee_surface.transform(geometry_context()).translation()(1),
422 perigee_surface.transform(geometry_context()).translation()(2));
423
424 trk.setChi2(gsftrk.chi2());
425 trk.setNhits(gsftrk.nMeasurements());
426 trk.setNdf(gsftrk.nMeasurements() - 5);
428 tracking::sim::utils::convertActsToLdmxPars(perigee_pars));
429 std::vector<double> v_trk_cov;
430 tracking::sim::utils::flatCov(trk_cov, v_trk_cov);
431 trk.setPerigeeCov(v_trk_cov);
432 Acts::Vector3 trk_momentum = gsftrk.momentum();
433 trk.setMomentum(trk_momentum(0), trk_momentum(1), trk_momentum(2));
434
435
436 trk.setTrackID(track.getTrackID());
437 trk.setPdgID(track.getPdgID());
438 trk.setTruthProb(track.getTruthProb());
439
440 itrk++;
441
442 out_tracks.push_back(trk);
443
444 }
445
446 event.add(out_trk_collection_, out_tracks);
447}
A source link that stores just an index.
bool exists(const std::string &name, const std::string &passName="", bool unique=true) const
Check for the existence of an object or collection with the given name and pass name in the event.
Implementation of a track object.
void setPerigeeParameters(const std::vector< double > &par)
d_0 z_0 phi_0 theta q/p t
void calibrate_1d(const Acts::GeometryContext &, const Acts::CalibrationContext &, const Acts::SourceLink &genericSourceLink, typename traj_t::TrackStateProxy trackState) const
Find the measurement corresponding to the source link.