Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Alignment.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Alignment.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2020-2021 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 
11 
12 template <typename fitter_t>
13 template <typename source_link_t, typename start_parameters_t,
14  typename fit_options_t>
18  const std::vector<source_link_t>& sourcelinks,
19  const start_parameters_t& sParameters, const fit_options_t& fitOptions,
20  const std::unordered_map<const Acts::Surface*, size_t>& idxedAlignSurfaces,
21  const ActsAlignment::AlignmentMask& alignMask) const {
24 
25  // Convert to Acts::SourceLink during iteration
26  Acts::SourceLinkAdapterIterator begin{sourcelinks.begin()};
27  Acts::SourceLinkAdapterIterator end{sourcelinks.end()};
28 
29  // Perform the fit
30  auto fitRes = m_fitter.fit(begin, end, sParameters, fitOptions, tracks);
31 
32  if (not fitRes.ok()) {
33  ACTS_WARNING("Fit failure");
34  return fitRes.error();
35  }
36  // The fit results
37  const auto& track = fitRes.value();
38  // Calculate the global track parameters covariance with the fitted track
39  const auto& globalTrackParamsCov =
41  tracks.trackStateContainer(), track.tipIndex());
42  // Calculate the alignment state
43  const auto alignState = detail::trackAlignmentState(
44  gctx, tracks.trackStateContainer(), track.tipIndex(),
45  globalTrackParamsCov, idxedAlignSurfaces, alignMask);
46  if (alignState.alignmentDof == 0) {
47  ACTS_VERBOSE("No alignment dof on track!");
48  return AlignmentError::NoAlignmentDofOnTrack;
49  }
50  return alignState;
51 }
52 
53 template <typename fitter_t>
54 template <typename trajectory_container_t,
55  typename start_parameters_container_t, typename fit_options_t>
57  const trajectory_container_t& trajectoryCollection,
58  const start_parameters_container_t& startParametersCollection,
59  const fit_options_t& fitOptions,
60  ActsAlignment::AlignmentResult& alignResult,
61  const ActsAlignment::AlignmentMask& alignMask) const {
62  // The number of trajectories must be equal to the number of starting
63  // parameters
64  assert(trajectoryCollection.size() == startParametersCollection.size());
65 
66  // The total alignment degree of freedom
67  alignResult.alignmentDof =
68  alignResult.idxedAlignSurfaces.size() * Acts::eAlignmentSize;
69  // Initialize derivative of chi2 w.r.t. alignment parameters for all tracks
70  Acts::ActsDynamicVector sumChi2Derivative =
71  Acts::ActsDynamicVector::Zero(alignResult.alignmentDof);
72  Acts::ActsDynamicMatrix sumChi2SecondDerivative =
73  Acts::ActsDynamicMatrix::Zero(alignResult.alignmentDof,
74  alignResult.alignmentDof);
75  // Copy the fit options
76  fit_options_t fitOptionsWithRefSurface = fitOptions;
77  // Calculate contribution to chi2 derivatives from all input trajectories
78  // @Todo: How to update the source link error iteratively?
79  alignResult.chi2 = 0;
80  alignResult.measurementDim = 0;
81  alignResult.numTracks = trajectoryCollection.size();
82  double sumChi2ONdf = 0;
83  for (unsigned int iTraj = 0; iTraj < trajectoryCollection.size(); iTraj++) {
84  const auto& sourcelinks = trajectoryCollection.at(iTraj);
85  const auto& sParameters = startParametersCollection.at(iTraj);
86  // Set the target surface
87  fitOptionsWithRefSurface.referenceSurface = &sParameters.referenceSurface();
88  // The result for one single track
89  auto evaluateRes = evaluateTrackAlignmentState(
90  fitOptions.geoContext, sourcelinks, sParameters,
91  fitOptionsWithRefSurface, alignResult.idxedAlignSurfaces, alignMask);
92  if (not evaluateRes.ok()) {
93  ACTS_DEBUG("Evaluation of alignment state for track " << iTraj
94  << " failed");
95  continue;
96  }
97  const auto& alignState = evaluateRes.value();
98  for (const auto& [rowSurface, rows] : alignState.alignedSurfaces) {
99  const auto& [dstRow, srcRow] = rows;
100  // Fill the results into full chi2 derivative matrix
101  sumChi2Derivative.segment<Acts::eAlignmentSize>(dstRow *
103  alignState.alignmentToChi2Derivative.segment(
104  srcRow * Acts::eAlignmentSize, Acts::eAlignmentSize);
105 
106  for (const auto& [colSurface, cols] : alignState.alignedSurfaces) {
107  const auto& [dstCol, srcCol] = cols;
108  sumChi2SecondDerivative
109  .block<Acts::eAlignmentSize, Acts::eAlignmentSize>(
110  dstRow * Acts::eAlignmentSize, dstCol * Acts::eAlignmentSize) +=
111  alignState.alignmentToChi2SecondDerivative.block(
112  srcRow * Acts::eAlignmentSize, srcCol * Acts::eAlignmentSize,
113  Acts::eAlignmentSize, Acts::eAlignmentSize);
114  }
115  }
116  alignResult.chi2 += alignState.chi2;
117  alignResult.measurementDim += alignState.measurementDim;
118  sumChi2ONdf += alignState.chi2 / alignState.measurementDim;
119  }
120  alignResult.averageChi2ONdf = sumChi2ONdf / alignResult.numTracks;
121 
122  // Get the inverse of chi2 second derivative matrix (we need this to
123  // calculate the covariance of the alignment parameters)
124  // @Todo: use more stable method for solving the inverse
125  size_t alignDof = alignResult.alignmentDof;
126  Acts::ActsDynamicMatrix sumChi2SecondDerivativeInverse =
127  Acts::ActsDynamicMatrix::Zero(alignDof, alignDof);
128  sumChi2SecondDerivativeInverse = sumChi2SecondDerivative.inverse();
129  if (sumChi2SecondDerivativeInverse.hasNaN()) {
130  ACTS_DEBUG("Chi2 second derivative inverse has NaN");
131  // return AlignmentError::AlignmentParametersUpdateFailure;
132  }
133 
134  // Initialize the alignment results
135  alignResult.deltaAlignmentParameters =
136  Acts::ActsDynamicVector::Zero(alignDof);
137  alignResult.alignmentCovariance =
138  Acts::ActsDynamicMatrix::Zero(alignDof, alignDof);
139  // Solve the linear equation to get alignment parameters change
140  alignResult.deltaAlignmentParameters =
141  -sumChi2SecondDerivative.fullPivLu().solve(sumChi2Derivative);
142  ACTS_VERBOSE("sumChi2SecondDerivative = \n" << sumChi2SecondDerivative);
143  ACTS_VERBOSE("sumChi2Derivative = \n" << sumChi2Derivative);
144  ACTS_VERBOSE("alignResult.deltaAlignmentParameters \n");
145 
146  // Alignment parameters covariance
147  alignResult.alignmentCovariance = 2 * sumChi2SecondDerivativeInverse;
148  // chi2 change
149  alignResult.deltaChi2 = 0.5 * sumChi2Derivative.transpose() *
150  alignResult.deltaAlignmentParameters;
151 }
152 
153 template <typename fitter_t>
156  const Acts::GeometryContext& gctx,
157  const std::vector<Acts::DetectorElementBase*>& alignedDetElements,
158  const ActsAlignment::AlignedTransformUpdater& alignedTransformUpdater,
159  ActsAlignment::AlignmentResult& alignResult) const {
160  // Update the aligned transform
161  Acts::AlignmentVector deltaAlignmentParam = Acts::AlignmentVector::Zero();
162  for (const auto& [surface, index] : alignResult.idxedAlignSurfaces) {
163  // 1. The original transform
164  const Acts::Vector3& oldCenter = surface->center(gctx);
165  const Acts::Transform3& oldTransform = surface->transform(gctx);
166  const Acts::RotationMatrix3& oldRotation = oldTransform.rotation();
167  // The elements stored below is (rotZ, rotY, rotX)
168  const Acts::Vector3& oldEulerAngles = oldRotation.eulerAngles(2, 1, 0);
169 
170  // 2. The delta transform
171  deltaAlignmentParam = alignResult.deltaAlignmentParameters.segment(
173  // The delta translation
174  Acts::Vector3 deltaCenter =
175  deltaAlignmentParam.segment<3>(Acts::eAlignmentCenter0);
176  // The delta Euler angles
177  Acts::Vector3 deltaEulerAngles =
178  deltaAlignmentParam.segment<3>(Acts::eAlignmentRotation0);
179 
180  // 3. The new transform
181  const Acts::Vector3 newCenter = oldCenter + deltaCenter;
182  // The rotation around global z axis
183  Acts::AngleAxis3 rotZ(oldEulerAngles(0) + deltaEulerAngles(2),
184  Acts::Vector3::UnitZ());
185  // The rotation around global y axis
186  Acts::AngleAxis3 rotY(oldEulerAngles(1) + deltaEulerAngles(1),
187  Acts::Vector3::UnitY());
188  // The rotation around global x axis
189  Acts::AngleAxis3 rotX(oldEulerAngles(2) + deltaEulerAngles(0),
190  Acts::Vector3::UnitX());
191  Eigen::Quaternion<Acts::ActsScalar> newRotation = rotZ * rotY * rotX;
192  const Acts::Transform3 newTransform =
193  Acts::Translation3(newCenter) * newRotation;
194 
195  // 4. Update the aligned transform
196  //@Todo: use a better way to handle this (need dynamic cast to inherited
197  // detector element type)
198  ACTS_VERBOSE("Delta of alignment parameters at element "
199  << index << "= \n"
200  << deltaAlignmentParam);
201  bool updated = alignedTransformUpdater(alignedDetElements.at(index), gctx,
202  newTransform);
203  if (not updated) {
204  ACTS_ERROR("Update alignment parameters for detector element failed");
205  return AlignmentError::AlignmentParametersUpdateFailure;
206  }
207  }
208 
210 }
211 
212 template <typename fitter_t>
213 template <typename trajectory_container_t,
214  typename start_parameters_container_t, typename fit_options_t>
217  const trajectory_container_t& trajectoryCollection,
218  const start_parameters_container_t& startParametersCollection,
219  const ActsAlignment::AlignmentOptions<fit_options_t>& alignOptions) const {
220  // Construct an AlignmentResult object
221  AlignmentResult alignResult;
222 
223  // Assign index to the alignable surface
224  for (unsigned int iDetElement = 0;
225  iDetElement < alignOptions.alignedDetElements.size(); iDetElement++) {
226  alignResult.idxedAlignSurfaces.emplace(
227  &alignOptions.alignedDetElements.at(iDetElement)->surface(),
228  iDetElement);
229  }
230  ACTS_VERBOSE("There are " << alignResult.idxedAlignSurfaces.size()
231  << " detector elements to be aligned");
232 
233  // Start the iteration to minimize the chi2
234  bool converged = false;
235  bool alignmentParametersUpdated = false;
236  std::queue<double> recentChi2ONdf;
237  ACTS_INFO("Max number of iterations: " << alignOptions.maxIterations);
238  for (unsigned int iIter = 0; iIter < alignOptions.maxIterations; iIter++) {
239  // Perform the fit to the trajectories and update alignment parameters
240  // Initialize the alignment mask (all dof in default)
241  AlignmentMask alignMask = AlignmentMask::All;
242  // Set the alignment mask
243  auto iter_it = alignOptions.iterationState.find(iIter);
244  if (iter_it != alignOptions.iterationState.end()) {
245  alignMask = iter_it->second;
246  }
247  // Calculate the alignment parameters delta etc.
248  calculateAlignmentParameters(
249  trajectoryCollection, startParametersCollection,
250  alignOptions.fitOptions, alignResult, alignMask);
251  // Screen out the information
252  ACTS_INFO("iIter = " << iIter << ", total chi2 = " << alignResult.chi2
253  << ", total measurementDim = "
254  << alignResult.measurementDim
255  << " and average chi2/ndf = "
256  << alignResult.averageChi2ONdf);
257  // Check if it has converged against the provided precision
258  // 1. either the delta average chi2/ndf in the last few
259  // iterations is within tolerance
260  if (recentChi2ONdf.size() >=
261  alignOptions.deltaAverageChi2ONdfCutOff.first) {
262  if (std::abs(recentChi2ONdf.front() - alignResult.averageChi2ONdf) <=
263  alignOptions.deltaAverageChi2ONdfCutOff.second) {
264  ACTS_INFO(
265  "Alignment has converged with change of chi2/ndf < "
266  << alignOptions.deltaAverageChi2ONdfCutOff.second << " in the last "
267  << alignOptions.deltaAverageChi2ONdfCutOff.first << " iterations"
268  << " after " << iIter << " iteration(s)");
269  converged = true;
270  break;
271  }
272  recentChi2ONdf.pop();
273  }
274  // 2. or the average chi2/ndf (is this correct?)
275  if (alignResult.averageChi2ONdf <= alignOptions.averageChi2ONdfCutOff) {
276  ACTS_INFO("Alignment has converged with average chi2/ndf < "
277  << alignOptions.averageChi2ONdfCutOff << " after " << iIter
278  << " iteration(s)");
279  converged = true;
280  break;
281  }
282  // Remove the first element
283  // Store the result in the queue
284  recentChi2ONdf.push(alignResult.averageChi2ONdf);
285 
286  ACTS_INFO("The solved delta of alignmentParameters = \n "
287  << alignResult.deltaAlignmentParameters);
288  // Not coveraged yet, update the detector element alignment parameters
289  auto updateRes = updateAlignmentParameters(
290  alignOptions.fitOptions.geoContext, alignOptions.alignedDetElements,
291  alignOptions.alignedTransformUpdater, alignResult);
292  if (not updateRes.ok()) {
293  ACTS_ERROR("Update alignment parameters failed: " << updateRes.error());
294  return updateRes.error();
295  }
296  alignmentParametersUpdated = true;
297  } // end of all iterations
298 
299  // Alignment failure if not converged
300  if (not converged) {
301  ACTS_ERROR("Alignment is not converged.");
302  alignResult.result = AlignmentError::ConvergeFailure;
303  }
304 
305  // Screen out the final aligned parameters
306  // @todo
307  if (alignmentParametersUpdated) {
308  for (const auto& det : alignOptions.alignedDetElements) {
309  const auto& surface = &det->surface();
310  const auto& transform =
311  det->transform(alignOptions.fitOptions.geoContext);
312  // write it to the result
313  alignResult.alignedParameters.emplace(det, transform);
314  const auto& translation = transform.translation();
315  const auto& rotation = transform.rotation();
316  const Acts::Vector3 rotAngles = rotation.eulerAngles(2, 1, 0);
317  ACTS_VERBOSE("Detector element with surface "
318  << surface->geometryId()
319  << " has aligned geometry position as below:");
320  ACTS_VERBOSE("Center (cenX, cenY, cenZ) = " << translation.transpose());
321  ACTS_VERBOSE(
322  "Euler angles (rotZ, rotY, rotX) = " << rotAngles.transpose());
323  ACTS_VERBOSE("Rotation matrix = \n" << rotation);
324  }
325  } else {
326  ACTS_DEBUG("Alignment parameters is not updated.");
327  }
328 
329  return alignResult;
330 }