TrackFittingAlgorithm.cpp 11.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
//
#include "TrackFittingAlgorithm.h"

// Gaudi
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/RndmGenerators.h"
10
#include "Gaudi/Property.h"
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"

#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
#include "Acts/Surfaces/PerigeeSurface.hpp"
#include "Acts/Propagator/EigenStepper.hpp"
#include "Acts/Propagator/Navigator.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Definitions/Common.hpp"
#include "Acts/Utilities/Helpers.hpp"
#include "Acts/Utilities/Logger.hpp"
#include "Acts/Definitions/Units.hpp"

#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
29
30
#include "JugBase/BField/DD4hepBField.h"

31
#include "JugTrack/GeometryContainers.hpp"
Whitney Armstrong's avatar
Whitney Armstrong committed
32
#include "JugTrack/IndexSourceLink.hpp"
33
34
#include "JugTrack/Track.hpp"
#include "JugTrack/Measurement.hpp"
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

#include "eicd/TrackerHitCollection.h"

#include <functional>
#include <stdexcept>
#include <vector>
#include <random>
#include <stdexcept>

namespace Jug::Reco {

  using namespace Acts::UnitLiterals;

  TrackFittingAlgorithm::TrackFittingAlgorithm(const std::string& name, ISvcLocator* svcLoc)
      : GaudiAlgorithm(name, svcLoc)
  {
    declareProperty("inputSourceLinks", m_inputSourceLinks, "");
    declareProperty("initialTrackParameters", m_initialTrackParameters, "");
Whitney Armstrong's avatar
Whitney Armstrong committed
53
    declareProperty("inputMeasurements", m_inputMeasurements, "");
54
    declareProperty("inputProtoTracks", m_inputProtoTracks, "");
55
56
57
58
59
60
61
62
63
64
65
66
67
68
    declareProperty("foundTracks", m_foundTracks, "");
    declareProperty("outputTrajectories", m_outputTrajectories, "");
  }

  StatusCode TrackFittingAlgorithm::initialize()
  {
    if (GaudiAlgorithm::initialize().isFailure())
      return StatusCode::FAILURE;
    m_geoSvc = service("GeoSvc");
    if (!m_geoSvc) {
      error() << "Unable to locate Geometry Service. "
              << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
      return StatusCode::FAILURE;
    }
69
70
    m_BField   = std::dynamic_pointer_cast<const Jug::BField::DD4hepBField>(m_geoSvc->getFieldProvider());
    m_fieldctx = Jug::BField::BFieldVariant(m_BField);
Whitney Armstrong's avatar
Whitney Armstrong committed
71

72
73
74
75
76
77
78
79
80
81
82
    // chi2 and #sourclinks per surface cutoffs
    //m_sourcelinkSelectorCfg = {
    //    {Acts::GeometryIdentifier(), {15, 10}},
    //};
    m_trackFittingFunc = makeTrackFittingFunction(m_geoSvc->trackingGeometry(), m_BField);
    return StatusCode::SUCCESS;
  }

  StatusCode TrackFittingAlgorithm::execute()
  {
    // Read input data
83
84
85
86
87
88
89
90
91
92
93
    const IndexSourceLinkContainer* sourceLinks       = m_inputSourceLinks.get();
    const TrackParametersContainer* initialParameters = m_initialTrackParameters.get();
    const MeasurementContainer*     measurements      = m_inputMeasurements.get();
    const ProtoTrackContainer*      protoTracks       = m_inputProtoTracks.get();
    ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("TrackFittingAlgorithm Logger", Acts::Logging::INFO));

    // Consistency cross checks
    if (protoTracks->size() != initialParameters->size()) {
      ACTS_FATAL("Inconsistent number of proto tracks and initial parameters");
      return StatusCode::FAILURE;
    }
94
95
96

    // TrajectoryContainer trajectories;
    auto trajectories = m_outputTrajectories.createAndPut();
97
    trajectories->reserve(initialParameters->size());
98

99
    // Construct a perigee surface as the target surface
100
101
    auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.});

102
103
    Acts::PropagatorPlainOptions pOptions;
    pOptions.maxSteps = 10000;
Whitney Armstrong's avatar
Whitney Armstrong committed
104
105
    // kfOptions.multipleScattering = m_cfg.multipleScattering;
    // kfOptions.energyLoss         = m_cfg.energyLoss;
106

107
    #if Acts_VERSION_MAJOR < 14
108
109
110
    Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder> kfOptions(
        m_geoctx, m_fieldctx, m_calibctx, MeasurementCalibrator(*measurements),
        Acts::VoidOutlierFinder(), Acts::LoggerWrapper{logger()}, pOptions, &(*pSurface));
111
112
113
114
115
116
    #else
    Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder, Acts::VoidReverseFilteringLogic> kfOptions(
        m_geoctx, m_fieldctx, m_calibctx, MeasurementCalibrator(*measurements),
        Acts::VoidOutlierFinder(), Acts::VoidReverseFilteringLogic(),
        Acts::LoggerWrapper{logger()}, pOptions, &(*pSurface));
    #endif
Whitney Armstrong's avatar
Whitney Armstrong committed
117

118
119
120
    // used for processing the data
    std::vector<IndexSourceLink>      trackSourceLinks;
    std::vector<const Acts::Surface*> surfSequence;
Whitney Armstrong's avatar
Whitney Armstrong committed
121

Sylvester Joosten's avatar
Sylvester Joosten committed
122
123
124
125
    if (msgLevel(MSG::DEBUG)) {
      debug() << "initialParams size:  " << initialParameters->size() << endmsg;
      debug() << "measurements size:  " << measurements->size() << endmsg;
    }
126

127
128
129
130
131
132
133
134
135
    // Perform the track finding for each starting parameter
    // @TODO: use seeds from track seeding algorithm as starting parameter
    // initial track params and proto tracks might likely have the same size.
    //for (std::size_t iseed = 0; iseed < init_trk_params->size(); ++iseed) {
    for (std::size_t itrack = 0; itrack < (*protoTracks).size(); ++itrack) {

      const auto& protoTrack    = (*protoTracks)[itrack];
      const auto& initialParams = (*initialParameters)[itrack];

Sylvester Joosten's avatar
Sylvester Joosten committed
136
137
138
139
      if (msgLevel(MSG::DEBUG)) {
        debug() << "protoTrack size:  " << protoTrack.size() << endmsg;
        debug() << "sourceLinks size:  " << sourceLinks->size() << endmsg;
      }
140
141
142
143
144

      trackSourceLinks.clear();
      trackSourceLinks.reserve(protoTrack.size());

      for (auto hitIndex : protoTrack) {
Sylvester Joosten's avatar
Sylvester Joosten committed
145
146
147
        if (msgLevel(MSG::DEBUG)) {
          debug() << " hit  index = " << hitIndex << endmsg;
        }
148
149
150
        auto sourceLink = sourceLinks->nth(hitIndex);
        auto geoId      = sourceLink->geometryId();
        if (sourceLink == sourceLinks->end()) {
Whitney Armstrong's avatar
Whitney Armstrong committed
151
          ACTS_FATAL("Proto track " << itrack << " contains invalid hit index "
152
153
154
155
156
157
                                    << hitIndex);
          return StatusCode::FAILURE;
        }
        trackSourceLinks.push_back(*sourceLink);
        //surfSequence.push_back(m_cfg.trackingGeometry->findSurface(geoId));
      }
158

Sylvester Joosten's avatar
Sylvester Joosten committed
159
160
161
      if (msgLevel(MSG::DEBUG)) {
        debug() << "Invoke track fitting ...  " << itrack << endmsg;
      }
Whitney Armstrong's avatar
Whitney Armstrong committed
162
      auto result = fitTrack(trackSourceLinks, initialParams, kfOptions);
Sylvester Joosten's avatar
Sylvester Joosten committed
163
164
165
      if (msgLevel(MSG::DEBUG)) {
        debug() << "fitting done." << endmsg;
      }
Whitney Armstrong's avatar
Whitney Armstrong committed
166
167
168
169
170
      // if (result.ok()) {
      //  // Get the track finding output object
      //  const auto& trackFindingOutput = result.value();
      //  // Create a SimMultiTrajectory
      //  trajectories->emplace_back(std::move(trackFindingOutput.fittedStates),
Wouter Deconinck's avatar
Acts v9    
Wouter Deconinck committed
171
      //  std::move(trackFindingOutput.lastMeasurementIndices),
Whitney Armstrong's avatar
Whitney Armstrong committed
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
      //                             std::move(trackFindingOutput.fittedParameters));
      //} else {
      //  debug() << "Track finding failed for truth seed " << iseed << endmsg;
      //  ACTS_WARNING("Track finding failed for truth seed " << iseed << " with error" <<
      //  result.error());
      //  // Track finding failed, but still create an empty SimMultiTrajectory
      //  // trajectories->push_back(SimMultiTrajectory());
      //}
      if (result.ok())
      {
        // Get the fit output object
        const auto& fitOutput = result.value();
        // The track entry indices container. One element here.
        std::vector<size_t> trackTips;
        trackTips.reserve(1);
        trackTips.emplace_back(fitOutput.lastMeasurementIndex);
        // The fitted parameters container. One element (at most) here.
        Trajectories::IndexedParameters indexedParams;
        //if (fitOutput.fittedParameters) {
        //  const auto& params = fitOutput.fittedParameters.value();
        //  ACTS_VERBOSE("Fitted paramemeters for track " << itrack);
        //  ACTS_VERBOSE("  " << params.parameters().transpose());
        //  // Push the fitted parameters to the container
        //  indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
        //} else {
        //  ACTS_DEBUG("No fitted paramemeters for track " << itrack);
        //}
        // store the result
        trajectories->emplace_back(std::move(fitOutput.fittedStates), std::move(trackTips),
                                  std::move(indexedParams));
202
      } else {
203
        ACTS_WARNING("Fit failed for track " << itrack << " with error" << result.error());
Whitney Armstrong's avatar
Whitney Armstrong committed
204
205
206
        // Fit failed. Add an empty result so the output container has
        // the same number of entries as the input.
        trajectories->push_back(Trajectories());
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
      }
    }

    // ctx.eventStore.add(m_cfg.outputTrajectories, std::move(trajectories));
    return StatusCode::SUCCESS;

    ///////////////////////////
    // acts example

  // Set the KalmanFitter options

  // Perform the fit for each input track
  //std::vector<IndexSourceLink> trackSourceLinks;
  //for (std::size_t itrack = 0; itrack < protoTracks.size(); ++itrack) {
  //  // The list of hits and the initial start parameters
  //  const auto& protoTrack = protoTracks[itrack];
  //  const auto& initialParams = initialParameters[itrack];

  //  // We can have empty tracks which must give empty fit results so the number
  //  // of entries in input and output containers matches.
  //  if (protoTrack.empty()) {
  //    trajectories.push_back(Trajectories());
  //    ACTS_WARNING("Empty track " << itrack << " found.");
  //    continue;
  //  }

  //  // Clear & reserve the right size
  //  trackSourceLinks.clear();
  //  trackSourceLinks.reserve(protoTrack.size());

  //  // Fill the source links via their indices from the container
  //  for (auto hitIndex : protoTrack) {
  //    auto sourceLink = sourceLinks.nth(hitIndex);
  //    if (sourceLink == sourceLinks.end()) {
  //      ACTS_FATAL("Proto track " << itrack << " contains invalid hit index"
  //                                << hitIndex);
  //      return ProcessCode::ABORT;
  //    }
  //    trackSourceLinks.push_back(*sourceLink);
  //  }

  //  ACTS_DEBUG("Invoke fitter");
  //  auto result = m_cfg.fit(trackSourceLinks, initialParams, kfOptions);
  //  if (result.ok()) {
  //    // Get the fit output object
  //    const auto& fitOutput = result.value();
  //    // The track entry indices container. One element here.
  //    std::vector<size_t> trackTips;
  //    trackTips.reserve(1);
  //    trackTips.emplace_back(fitOutput.trackTip);
  //    // The fitted parameters container. One element (at most) here.
  //    Trajectories::IndexedParameters indexedParams;
  //    if (fitOutput.fittedParameters) {
  //      const auto& params = fitOutput.fittedParameters.value();
  //      ACTS_VERBOSE("Fitted paramemeters for track " << itrack);
  //      ACTS_VERBOSE("  " << params.parameters().transpose());
  //      // Push the fitted parameters to the container
  //      indexedParams.emplace(fitOutput.trackTip, std::move(params));
  //    } else {
  //      ACTS_DEBUG("No fitted paramemeters for track " << itrack);
  //    }
  //    // store the result
  //    trajectories.emplace_back(std::move(fitOutput.fittedStates),
  //                              std::move(trackTips), std::move(indexedParams));
  //  } else {
  //    ACTS_WARNING("Fit failed for track " << itrack << " with error"
  //                                         << result.error());
  //    // Fit failed. Add an empty result so the output container has
  //    // the same number of entries as the input.
  //    trajectories.push_back(Trajectories());
  //  }
  //}


    return StatusCode::SUCCESS;
  }

  DECLARE_COMPONENT(TrackFittingAlgorithm)
} // namespace Jug::Reco