IdealResult Class Reference

Calculates the ideal result of the StraightTrackAlignment. More...

#include <IdealResult.h>

Inheritance diagram for IdealResult:

Public Member Functions
void	analyze () override
	analyzes the data collected More...
void	begin (const CTPPSGeometry *geometryReal, const CTPPSGeometry *geometryMisaligned) override
	prepare for processing More...
void	end () override
	cleans up after processing More...
void	feed (const HitCollection &, const LocalTrackFit &) override
	process one track More...
std::string	getName () override
bool	hasErrorEstimate () override
	returns whether this algorithm is capable of estimating result uncertainties More...
	IdealResult ()
	dummy constructor (not to be used) More...
	IdealResult (const edm::ParameterSet &ps, AlignmentTask *_t)
	normal constructor More...
void	saveDiagnostics (TDirectory *) override
	saves diagnostic histograms/plots More...
unsigned int	solve (const std::vector< AlignmentConstraint > &, std::map< unsigned int, AlignmentResult > &result, TDirectory *dir) override
	~IdealResult () override
Public Member Functions inherited from AlignmentAlgorithm
	AlignmentAlgorithm ()
	dummy constructor (not to be used) More...
	AlignmentAlgorithm (const edm::ParameterSet &ps, AlignmentTask *_t)
	normal constructor More...
virtual	~AlignmentAlgorithm ()

Protected Attributes
edm::ESHandle< CTPPSGeometry >	gMisaligned
edm::ESHandle< CTPPSGeometry >	gReal
Protected Attributes inherited from AlignmentAlgorithm
double	singularLimit
	eigenvalues in (-singularLimit, singularLimit) are treated as singular More...
AlignmentTask *	task
	the tasked to be completed More...
unsigned int	verbosity

Detailed Description

Calculates the ideal result of the StraightTrackAlignment.

Definition at line 24 of file IdealResult.h.

Constructor & Destructor Documentation

IdealResult::IdealResult ( )

inline

dummy constructor (not to be used)

Definition at line 30 of file IdealResult.h.

{}

IdealResult::IdealResult	(	const edm::ParameterSet &	ps,
		AlignmentTask *	_t
	)

normal constructor

Definition at line 20 of file IdealResult.cc.

: AlignmentAlgorithm(gps, _t) {}

IdealResult::~IdealResult ( )

inlineoverride

Definition at line 35 of file IdealResult.h.

{}

Member Function Documentation

void IdealResult::analyze ( )

inlineoverridevirtual

analyzes the data collected

Implements AlignmentAlgorithm.

Definition at line 47 of file IdealResult.h.

{}

void IdealResult::begin ( const CTPPSGeometry *  geometryReal, const CTPPSGeometry *  geometryMisaligned  )

overridevirtual

prepare for processing

Implements AlignmentAlgorithm.

Definition at line 24 of file IdealResult.cc.

References gMisaligned, and gReal.

                                                                                                  {
  gReal = geometryReal;
  gMisaligned = geometryMisaligned;
}

void IdealResult::end ( )

inlineoverridevirtual

cleans up after processing

Implements AlignmentAlgorithm.

Definition at line 53 of file IdealResult.h.

Referenced by Types.LuminosityBlockRange::cppID(), and Types.EventRange::cppID().

{}

void IdealResult::feed ( const HitCollection &  , const LocalTrackFit &    )

inlineoverridevirtual

process one track

Implements AlignmentAlgorithm.

Definition at line 43 of file IdealResult.h.

{}

std::string IdealResult::getName ( )

inlineoverridevirtual

Reimplemented from AlignmentAlgorithm.

Definition at line 37 of file IdealResult.h.

Referenced by plotting.Plot::draw().

{ return "Ideal"; }

bool IdealResult::hasErrorEstimate ( )

inlineoverridevirtual

returns whether this algorithm is capable of estimating result uncertainties

Implements AlignmentAlgorithm.

Definition at line 39 of file IdealResult.h.

{ return false; }

void IdealResult::saveDiagnostics ( TDirectory *  )

inlineoverridevirtual

saves diagnostic histograms/plots

Implements AlignmentAlgorithm.

Definition at line 45 of file IdealResult.h.

{}

unsigned int IdealResult::solve ( const std::vector< AlignmentConstraint > &  , std::map< unsigned int, AlignmentResult > &  results, TDirectory *  dir  )

overridevirtual

solves the alignment problem with the given constraints

Parameters

dir	a directory (in StraightTrackAlignment::taskDataFileName) where intermediate results can be stored

Implements AlignmentAlgorithm.

Definition at line 30 of file IdealResult.cc.

References funct::A, funct::abs(), b, gen::C, cuy::cv, ztail::d, Exception, relativeConstraints::geom, AlignmentTask::geometry, AlignmentGeometry::get(), AlignmentTask::getQuantityIndex(), AlignmentGeometry::getSensorMap(), gMisaligned, gReal, Exhume::I, mps_fire::i, dqmiolumiharvest::j, isotrackApplyRegressor::k, visualization-live-secondInstance_cfg::m, makeMuonMisalignmentScenario::misal, hltrates_dqm_sourceclient-live_cfg::offset, upgradeWorkflowComponents::offsets, AlignmentTask::qcRotZ, AlignmentTask::qcShR1, AlignmentTask::qcShR2, AlignmentTask::qcShZ, AlignmentTask::quantitiesOfClass(), AlignmentTask::quantityClasses, alignCSCRings::r, AlignmentTask::resolveRotZ, AlignmentTask::resolveShR, rho, idealTransformation::rotation, DetGeomDesc::rotation(), CTPPSDetId::sdTrackingPixel, CTPPSDetId::sdTrackingStrip, AlignmentResult::setRotZ(), AlignmentResult::setShR1(), AlignmentResult::setShR2(), AlignmentResult::setShZ(), edm::shift, DetId::subdetId(), AlignmentAlgorithm::task, DetGeomDesc::translation(), cms::cuda::V, and findQualityFiles::v.

                                                 {
  /* 
  STRATEGY:

  1) Determine true misalignment as misalign_geometry - real_geometry.

  2) Represent the true misalignment as the "chi" vector (cf. Jan algorithm), denoted chi_tr.

  3) Find the expected result of track-based alignment, subject to the given constraints. Formally this corresponds to finding
       chi_exp = chi_tr + I * Lambda
     where the I matrix contains the "inaccessible alignment modes" as columns and Lambda is a vector of parameters. The constraints are given by
       C^T * chi_exp = V
     Since the problem may be generally overconstrained, it is better to formulate it as search for Lambda which minimises
       |C^ * chi_exp - V|^2
     This gives
       Lambda = (A^T * A)^-1 * A^T * b,  A = C^T * I,  b = V - C^T * chi_tr
  */

  results.clear();

  // determine dimension and offsets
  unsigned int dim = 0;
  vector<unsigned int> offsets;
  map<unsigned int, unsigned int> offset_map;
  for (unsigned int qci = 0; qci < task->quantityClasses.size(); qci++) {
    offsets.push_back(dim);
    offset_map[task->quantityClasses[qci]] = dim;
    dim += task->quantitiesOfClass(task->quantityClasses[qci]);
  }

  // collect true misalignments
  TVectorD chi_tr(dim);

  for (const auto &dit : task->geometry.getSensorMap()) {
    unsigned int detId = dit.first;

    const DetGeomDesc *real = gReal->sensor(detId);
    const DetGeomDesc *misal = gMisaligned->sensor(detId);

    // extract shift
    const auto shift = misal->translation() - real->translation();

    // extract rotation around z
    const auto rotation = misal->rotation() * real->rotation().Inverse();

    double r_xx, r_xy, r_xz;
    double r_yx, r_yy, r_yz;
    double r_zx, r_zy, r_zz;
    rotation.GetComponents(r_xx, r_xy, r_xz, r_yx, r_yy, r_yz, r_zx, r_zy, r_zz);

    if (std::abs(r_zz - 1.) > 1E-5)
      throw cms::Exception("PPS") << "IdealResult::Solve: only rotations about z are supported.";

    double rot_z = atan2(r_yx, r_xx);

    const auto &geom = task->geometry.get(detId);

    for (unsigned int qci = 0; qci < task->quantityClasses.size(); ++qci) {
      const auto &qc = task->quantityClasses[qci];
      signed int idx = task->getQuantityIndex(qc, detId);
      if (idx < 0)
        continue;

      double v = 0.;

      if (qc == AlignmentTask::qcShR1) {
        const auto &d = geom.getDirectionData(1);
        v = shift.x() * d.dx + shift.y() * d.dy + shift.z() * d.dz;
      }

      if (qc == AlignmentTask::qcShR2) {
        const auto &d = geom.getDirectionData(2);
        v = shift.x() * d.dx + shift.y() * d.dy + shift.z() * d.dz;
      }

      if (qc == AlignmentTask::qcRotZ)
        v = rot_z;

      chi_tr(offsets[qci] + idx) = v;
    }
  }

  // build list of "inaccessible" modes
  vector<TVectorD> inaccessibleModes;

  if (task->resolveShR) {
    TVectorD fm_ShX_gl(dim);
    fm_ShX_gl.Zero();
    TVectorD fm_ShX_lp(dim);
    fm_ShX_lp.Zero();
    TVectorD fm_ShY_gl(dim);
    fm_ShY_gl.Zero();
    TVectorD fm_ShY_lp(dim);
    fm_ShY_lp.Zero();

    for (const auto &sp : task->geometry.getSensorMap()) {
      CTPPSDetId senId(sp.first);
      const auto &geom = sp.second;

      if (senId.subdetId() == CTPPSDetId::sdTrackingStrip) {
        signed int qIndex = task->getQuantityIndex(AlignmentTask::qcShR2, senId);

        const double d2x = geom.getDirectionData(2).dx;
        const double d2y = geom.getDirectionData(2).dy;

        const auto &offset2 = offset_map[AlignmentTask::qcShR2];
        fm_ShX_gl(offset2 + qIndex) = d2x;
        fm_ShX_lp(offset2 + qIndex) = d2x * geom.z;
        fm_ShY_gl(offset2 + qIndex) = d2y;
        fm_ShY_lp(offset2 + qIndex) = d2y * geom.z;
      }

      if (senId.subdetId() == CTPPSDetId::sdTrackingPixel) {
        const signed int qIndex1 = task->getQuantityIndex(AlignmentTask::qcShR1, senId);
        const signed int qIndex2 = task->getQuantityIndex(AlignmentTask::qcShR2, senId);

        const double d1x = geom.getDirectionData(1).dx;
        const double d1y = geom.getDirectionData(1).dy;
        const double d2x = geom.getDirectionData(2).dx;
        const double d2y = geom.getDirectionData(2).dy;

        const auto &offset1 = offset_map[AlignmentTask::qcShR1];
        fm_ShX_gl(offset1 + qIndex1) = d1x;
        fm_ShX_lp(offset1 + qIndex1) = d1x * geom.z;
        fm_ShY_gl(offset1 + qIndex1) = d1y;
        fm_ShY_lp(offset1 + qIndex1) = d1y * geom.z;

        const auto &offset2 = offset_map[AlignmentTask::qcShR2];
        fm_ShX_gl(offset2 + qIndex2) = d2x;
        fm_ShX_lp(offset2 + qIndex2) = d2x * geom.z;
        fm_ShY_gl(offset2 + qIndex2) = d2y;
        fm_ShY_lp(offset2 + qIndex2) = d2y * geom.z;
      }
    }

    inaccessibleModes.push_back(fm_ShX_gl);
    inaccessibleModes.push_back(fm_ShX_lp);
    inaccessibleModes.push_back(fm_ShY_gl);
    inaccessibleModes.push_back(fm_ShY_lp);
  }

  if (task->resolveRotZ) {
    TVectorD fm_RotZ_gl(dim);
    fm_RotZ_gl.Zero();
    TVectorD fm_RotZ_lp(dim);
    fm_RotZ_lp.Zero();

    for (const auto &sp : task->geometry.getSensorMap()) {
      CTPPSDetId senId(sp.first);
      const auto &geom = sp.second;

      for (int m = 0; m < 2; ++m) {
        double rho = 0.;
        TVectorD *fmp = nullptr;
        if (m == 0) {
          rho = 1.;
          fmp = &fm_RotZ_gl;
        }
        if (m == 1) {
          rho = geom.z;
          fmp = &fm_RotZ_lp;
        }
        TVectorD &fm = *fmp;

        const signed int qIndex = task->getQuantityIndex(AlignmentTask::qcRotZ, senId);
        const auto &offset = offset_map[AlignmentTask::qcRotZ];
        fm(offset + qIndex) = rho;

        const double as_x = -rho * geom.sy;
        const double as_y = +rho * geom.sx;

        if (senId.subdetId() == CTPPSDetId::sdTrackingStrip) {
          const double d2x = geom.getDirectionData(2).dx;
          const double d2y = geom.getDirectionData(2).dy;

          const signed int qIndex2 = task->getQuantityIndex(AlignmentTask::qcShR2, senId);
          const auto &offset2 = offset_map[AlignmentTask::qcShR2];
          fm(offset2 + qIndex2) = d2x * as_x + d2y * as_y;
        }

        if (senId.subdetId() == CTPPSDetId::sdTrackingPixel) {
          const double d1x = geom.getDirectionData(1).dx;
          const double d1y = geom.getDirectionData(1).dy;
          const double d2x = geom.getDirectionData(2).dx;
          const double d2y = geom.getDirectionData(2).dy;

          const signed int qIndex1 = task->getQuantityIndex(AlignmentTask::qcShR1, senId);
          const auto &offset1 = offset_map[AlignmentTask::qcShR1];
          fm(offset1 + qIndex1) = d1x * as_x + d1y * as_y;

          const signed int qIndex2 = task->getQuantityIndex(AlignmentTask::qcShR2, senId);
          const auto &offset2 = offset_map[AlignmentTask::qcShR2];
          fm(offset2 + qIndex2) = d2x * as_x + d2y * as_y;
        }
      }
    }

    inaccessibleModes.push_back(fm_RotZ_gl);
    inaccessibleModes.push_back(fm_RotZ_lp);
  }

  // build matrices and vectors
  TMatrixD C(dim, constraints.size());
  TVectorD V(constraints.size());
  for (unsigned int i = 0; i < constraints.size(); i++) {
    V(i) = constraints[i].val;

    unsigned int offset = 0;
    for (auto &quantityClass : task->quantityClasses) {
      const TVectorD &cv = constraints[i].coef.find(quantityClass)->second;
      for (int k = 0; k < cv.GetNrows(); k++) {
        C[offset][i] = cv[k];
        offset++;
      }
    }
  }

  TMatrixD I(dim, inaccessibleModes.size());
  for (unsigned int i = 0; i < inaccessibleModes.size(); ++i) {
    for (int j = 0; j < inaccessibleModes[i].GetNrows(); ++j)
      I(j, i) = inaccessibleModes[i](j);
  }

  // determine expected track-based alignment result
  TMatrixD CT(TMatrixD::kTransposed, C);
  TMatrixD CTI(CT * I);

  const TMatrixD &A = CTI;
  TMatrixD AT(TMatrixD::kTransposed, A);
  TMatrixD ATA(AT * A);
  TMatrixD ATA_inv(TMatrixD::kInverted, ATA);

  TVectorD b = V - CT * chi_tr;

  TVectorD La(ATA_inv * AT * b);

  TVectorD chi_exp(chi_tr + I * La);

  // save result
  for (const auto &dit : task->geometry.getSensorMap()) {
    AlignmentResult r;

    for (unsigned int qci = 0; qci < task->quantityClasses.size(); ++qci) {
      const auto &qc = task->quantityClasses[qci];
      const auto idx = task->getQuantityIndex(qc, dit.first);
      if (idx < 0)
        continue;

      const auto &v = chi_exp(offsets[qci] + idx);

      switch (qc) {
        case AlignmentTask::qcShR1:
          r.setShR1(v);
          break;
        case AlignmentTask::qcShR2:
          r.setShR2(v);
          break;
        case AlignmentTask::qcShZ:
          r.setShZ(v);
          break;
        case AlignmentTask::qcRotZ:
          r.setRotZ(v);
          break;
      }
    }

    results[dit.first] = r;
  }

  return 0;
}

Member Data Documentation

edm::ESHandle<CTPPSGeometry> IdealResult::gMisaligned

protected

Definition at line 26 of file IdealResult.h.

Referenced by begin(), and solve().

edm::ESHandle<CTPPSGeometry> IdealResult::gReal

protected

Definition at line 26 of file IdealResult.h.

Referenced by begin(), and solve().