Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Alignment/CommonAlignment/interface/Alignable.h"
#include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
#include "Alignment/CommonAlignment/interface/AlignmentParameters.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"

#include "Alignment/CommonAlignmentParametrization/interface/CompositeAlignmentDerivativesExtractor.h"

//--------------------------------------------------------------------------------------
CompositeAlignmentDerivativesExtractor::CompositeAlignmentDerivativesExtractor(
    const align::Alignables &alignables,
    const std::vector<AlignableDet *> &alignableDets,
    const std::vector<TrajectoryStateOnSurface> &tsos) {
  std::vector<AlignableDetOrUnitPtr> detOrUnits;
  detOrUnits.reserve(alignableDets.size());

  std::vector<AlignableDet *>::const_iterator it, itEnd;
  for (it = alignableDets.begin(), itEnd = alignableDets.end(); it != itEnd; ++it)
    detOrUnits.push_back(AlignableDetOrUnitPtr(*it));

  extractCurrentAlignment(alignables, detOrUnits, tsos);
}

//--------------------------------------------------------------------------------------
CompositeAlignmentDerivativesExtractor::CompositeAlignmentDerivativesExtractor(
    const align::Alignables &alignables,
    const std::vector<AlignableDetOrUnitPtr> &alignableDets,
    const std::vector<TrajectoryStateOnSurface> &tsos) {
  extractCurrentAlignment(alignables, alignableDets, tsos);
}

//--------------------------------------------------------------------------------------

void CompositeAlignmentDerivativesExtractor::extractCurrentAlignment(
    const align::Alignables &alignables,
    const std::vector<AlignableDetOrUnitPtr> &alignableDets,
    const std::vector<TrajectoryStateOnSurface> &tsos) {
  // sanity check
  if (alignables.size() != alignableDets.size()) {
    edm::LogError("CompositeAlignmentDerivativesExtractor")
        << "Inconsistent length of arguments: alignables=" << alignables.size()
        << ", alignableDets=" << alignableDets.size();
    return;
  }

  if (alignables.size() != tsos.size()) {
    edm::LogError("CompositeAlignmentDerivativesExtractor")
        << "Inconsistent length of arguments: alignables=" << alignables.size() << ", tsos=" << tsos.size();
    return;
  }

  align::Alignables::const_iterator itAlignable = alignables.begin();
  std::vector<AlignableDetOrUnitPtr>::const_iterator itAlignableDet = alignableDets.begin();
  std::vector<TrajectoryStateOnSurface>::const_iterator itTsos = tsos.begin();

  int nRow = 0;
  int nCollumn = 0;
  unsigned int nAlignables = 0;

  std::vector<AlgebraicMatrix> subDerivatives;
  std::vector<AlgebraicVector> subCorrectionTerm;

  // get the individual derivatives and correction term and determine the
  // dimension
  while (itAlignable != alignables.end()) {
    // Get the current estimate on the alignment parameters
    AlgebraicVector subAlignmentParameters = (*itAlignable)->alignmentParameters()->selectedParameters();

    // Get the derivatives or the local coordinates w.r.t. the corresponding
    // alignment parameters
    AlgebraicMatrix subAlignmentDerivatives =
        (*itAlignable)->alignmentParameters()->selectedDerivatives(*itTsos, *itAlignableDet);

    subDerivatives.push_back(subAlignmentDerivatives.T());
    subCorrectionTerm.push_back(subAlignmentDerivatives.T() * subAlignmentParameters);

    nRow += 2;
    // check if it is the first occurrence of this Alignable
    if (count(alignables.begin(), itAlignable, *itAlignable) == 0) {
      // matrix is transposed -> num_row() instead of num_col()
      nCollumn += subAlignmentDerivatives.num_row();
      nAlignables++;
    }

    ++itAlignable;
    ++itAlignableDet;
    ++itTsos;
  }

  // construct derivatives and correction term with the right dimension
  theDerivatives = AlgebraicMatrix(nRow, nCollumn, 0);
  theCorrectionTerm = AlgebraicVector(nRow, 0);

  if (alignables.size() == nAlignables)
    // One hit per alignable
    extractWithoutMultipleHits(subCorrectionTerm, subDerivatives);
  else
    // At least one alignable has two hits
    extractWithMultipleHits(subCorrectionTerm, subDerivatives, alignables);

  return;
}

//--------------------------------------------------------------------------------------

void CompositeAlignmentDerivativesExtractor::extractWithoutMultipleHits(
    const std::vector<AlgebraicVector> &subCorrectionTerm, const std::vector<AlgebraicMatrix> &subDerivatives) {
  std::vector<AlgebraicVector>::const_iterator itSubCorrectionTerm = subCorrectionTerm.begin();
  std::vector<AlgebraicMatrix>::const_iterator itSubDerivatives = subDerivatives.begin();

  int iRow = 1;
  int iCollumn = 1;

  // Fill in the individual terms
  while (itSubCorrectionTerm != subCorrectionTerm.end()) {
    theCorrectionTerm.sub(iRow, *itSubCorrectionTerm);
    theDerivatives.sub(iRow, iCollumn, *itSubDerivatives);

    iRow += 2;
    iCollumn += (*itSubDerivatives).num_col();

    ++itSubCorrectionTerm;
    ++itSubDerivatives;
  }

  return;
}

//--------------------------------------------------------------------------------------

void CompositeAlignmentDerivativesExtractor::extractWithMultipleHits(
    const std::vector<AlgebraicVector> &subCorrectionTerm,
    const std::vector<AlgebraicMatrix> &subDerivatives,
    const align::Alignables &alignables) {
  std::vector<AlgebraicVector>::const_iterator itSubCorrectionTerm = subCorrectionTerm.begin();
  std::vector<AlgebraicMatrix>::const_iterator itSubDerivatives = subDerivatives.begin();
  align::Alignables::const_iterator itAlignables = alignables.begin();
  align::Alignables::const_iterator itPosition;
  align::Alignables::const_iterator itLastPosition;

  int iRow = 1;

  // Fill in the individual terms
  while (itAlignables != alignables.end()) {
    theCorrectionTerm.sub(iRow, *itSubCorrectionTerm);

    int iCollumn = 1;
    int iAlignable = 0;

    itLastPosition = find(alignables.begin(), itAlignables, *itAlignables);

    for (itPosition = alignables.begin(); itPosition != itLastPosition; ++itPosition) {
      if (count(alignables.begin(), itPosition, *itPosition) == 0)
        iCollumn += subDerivatives[iAlignable].num_col();
      iAlignable++;
    }

    theDerivatives.sub(iRow, iCollumn, *itSubDerivatives);

    iRow += 2;

    ++itAlignables;
    ++itSubCorrectionTerm;
    ++itSubDerivatives;
  }

  return;
}