PPS: relative alignment (#31629)

* Added PPS relative alignment.

* Fixed logic in getFullSensorCorrection.

* Apply code check and format patches.

* Changes suggested by code checks.

* Added invertibility checks.

* Removed obsolete code.

* Use module types from threaded framework.

* Avoid memory-leak warning.

* Avoid logic-error warning.

* Use es token, part one.

* Use es token, part two.

* Clang format.

* Catch only cms::Exception.

CalibPPS/AlignmentRelative/BuildFile.xml

@@ -0,0 +1,13 @@
+<use name="clhep"/>
+<use name="root"/>
+<use name="rootgraphics"/>
+
+<use name="DataFormats/CTPPSDetId"/>
+<use name="Geometry/VeryForwardGeometryBuilder"/>
+<use name="Geometry/Records"/>
+<use name="DataFormats/CTPPSReco"/>
+<use name="CondFormats/PPSObjects"/>
+
+<export>
+ <lib name="1"/>
+</export>

CalibPPS/AlignmentRelative/README.md

@@ -0,0 +1,61 @@
+# introduction
+
+Alignment of (CT-)PPS detectors (housed in Roman Pots) proceeds in the following conceptual steps
+ * relative alignment among sensors - by minimisation of track-hit residuals
+ * global alignment of RPs wrt the beam - by exploiting symmetries in observed hit patterns
+
+The alignment is applied separately to each arm of the spectrometer (LHC sectors 45 adn 56). For more details see e.g. this [note](http://cds.cern.ch/record/2256296).
+
+This package implements the first step, "track-based alignment" among RP sensors. The alignment corrections of relevance are transverse shifts (in x and y) a rotations about the beam axis (z). The method is based on minimisation of residuals between hits and reconstructed tracks. Certain misalignment modes do not generate residuals (e.g. global shift or rotation) and are therefore inaccessible to the track-based alignment. Sometimes, these modes are also referred to as "singular". In order to find a solution to the alignment task, user must specify "constraints" which provide information about the inaccessible/singular alignment modes. These constraints may come from the second step (global alignment) as outlined above.
+
+For the following reasons, the track-based alignment is performed in an iterative manner:
+ * the treatment of rotations in linearised (sin alpha ~ alpha, etc.)
+ * a priory, it is not clear that a hit with large residual is either an outlier (noise or so) or hit from a sensor with large misalignment
+
+Therefore the alignment starts with a large tolerace (don't exclude hits with large misalignments) and gradually (in iterations over the same data) makes the quality cuts stricter (to remove outliers).
+
+The implementation inherits from the code originally developed for TOTEM and described in this [thesis](http://cdsweb.cern.ch/record/1441140). The original code was developed for Si strip detectors, the current implementation can cope also with the Si pixel detectors of PPS. There is some code also for the timing RPs (diamonds), but it is of experimental nature.
+
+# input
+
+For strip RPs, the alignment takes the input from "U-V patterns", i.e. linear rec-hit patterns in U-z and V-z planes.
+
+For pixel RPs, the alignment can take the input
+ * either from rec-hits (only from planes with a single hit)
+ * or reconstructed tracks (reconstruction provides some suppression of outliers).
+
+One should not use both input methods in the same time (double counting). 
+
+# files
+
+In interface/
+ * AlignmentAlgorithm.h: abstract interface of a track-based algorithm
+ * AlignmentConstraint.h: information about a constraint (fixing an alignment mode inaccessible to track-based alignment)
+ * AlignmentGeometry.h: summary of geometry-related data used in alignment analysis
+ * AlignmentResult.h: structure holding the track-based alignment results
+ * AlignmentTask.h: information about alignment task
+ * HitCollection.h: collection of sensor hits (from strip or pixel sensors)
+ * IdealResult.h: predicts the result of track-based alignment under the imposed constraints
+ * JanAlignmentAlgorithm.h: an implementation of track-hit minimisation algorithm
+ * LocalTrackFit.h: straight-line fit through all RPs in one arm
+ * LocalTrackFitter.h: code to make the per-arm track fits
+ * SingularMode.h: information about inaccessible/singular modes
+ * StraightTrackAlignment.h: common code for all track-based alignment algorithms
+ * Utilities.h: common code
+
+In plugins/
+ * PPSFastLocalSimulation.cc: a fast per-arm simulation to test the alignment code
+ * PPSModifySingularModes.cc: module to modify the singular modes in a track-based alignment result
+ * PPSStraightTrackAligner.cc: module to run the track-based alignment
+
+# tests
+
+ * test_with_mc: Monte-Carlo tests of the alignment code
+   * simple: a very basic test
+   * iterations: script to test alignment convergence in several iterations
+   * statistics: script to test statistical properties
+
+ * test_modify_singular_modes: an example how to modify the singular modes
+
+ * test_with_data: a full/real-life example of alignment application to (LHC) data
+

CalibPPS/AlignmentRelative/interface/AlignmentAlgorithm.h

@@ -0,0 +1,78 @@
+/****************************************************************************
+* Authors: 
+*  Jan Kašpar ([email protected]) 
+****************************************************************************/
+
+#ifndef CalibPPS_AlignmentRelative_AlignmentAlgorithm_h
+#define CalibPPS_AlignmentRelative_AlignmentAlgorithm_h
+
+#include "FWCore/Framework/interface/EventSetup.h"
+
+#include "CalibPPS/AlignmentRelative/interface/LocalTrackFit.h"
+#include "CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h"
+#include "CalibPPS/AlignmentRelative/interface/HitCollection.h"
+#include "CalibPPS/AlignmentRelative/interface/AlignmentConstraint.h"
+#include "CalibPPS/AlignmentRelative/interface/AlignmentResult.h"
+#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
+
+#include <string>
+#include <map>
+
+class AlignmentTask;
+class TDirectory;
+
+namespace edm {
+  class ParameterSet;
+}
+
+/**
+ *\brief Abstract parent for all (track-based) alignment algorithms
+ **/
+class AlignmentAlgorithm {
+protected:
+  unsigned int verbosity;
+
+  /// the tasked to be completed
+  AlignmentTask *task;
+
+  /// eigenvalues in (-singularLimit, singularLimit) are treated as singular
+  double singularLimit;
+
+public:
+  /// dummy constructor (not to be used)
+  AlignmentAlgorithm() {}
+
+  /// normal constructor
+  AlignmentAlgorithm(const edm::ParameterSet &ps, AlignmentTask *_t);
+
+  virtual ~AlignmentAlgorithm() {}
+
+  virtual std::string getName() { return "Base"; }
+
+  /// returns whether this algorithm is capable of estimating result uncertainties
+  virtual bool hasErrorEstimate() = 0;
+
+  /// prepare for processing
+  virtual void begin(const CTPPSGeometry *geometryReal, const CTPPSGeometry *geometryMisaligned) = 0;
+
+  /// process one track
+  virtual void feed(const HitCollection &, const LocalTrackFit &) = 0;
+
+  /// saves diagnostic histograms/plots
+  virtual void saveDiagnostics(TDirectory *) = 0;
+
+  /// analyzes the data collected
+  virtual void analyze() = 0;
+
+  /// solves the alignment problem with the given constraints
+  /// \param dir a directory (in StraightTrackAlignment::taskDataFileName) where
+  /// intermediate results can be stored
+  virtual unsigned int solve(const std::vector<AlignmentConstraint> &,
+                             std::map<unsigned int, AlignmentResult> &results,
+                             TDirectory *dir = nullptr) = 0;
+
+  /// cleans up after processing
+  virtual void end() = 0;
+};
+
+#endif

CalibPPS/AlignmentRelative/interface/AlignmentConstraint.h

@@ -0,0 +1,29 @@
+/****************************************************************************
+* Authors: 
+*  Jan Kašpar ([email protected]) 
+****************************************************************************/
+
+#ifndef CalibPPS_AlignmentRelative_AlignmentConstraint_h
+#define CalibPPS_AlignmentRelative_AlignmentConstraint_h
+
+#include <TVectorD.h>
+
+#include <map>
+#include <string>
+
+/**
+ *\brief An alignment constraint.
+ **/
+class AlignmentConstraint {
+public:
+  /// constraint value
+  double val;
+
+  /// map: AlignmentAlgorithm::QuantityClass -> constraint coefficients
+  std::map<unsigned int, TVectorD> coef;
+
+  /// label of the constraint
+  std::string name;
+};
+
+#endif

CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h

@@ -0,0 +1,82 @@
+/****************************************************************************
+* Authors: 
+*  Jan Kašpar ([email protected]) 
+****************************************************************************/
+
+#ifndef CalibPPS_AlignmentRelative_AlignmentGeometry_h
+#define CalibPPS_AlignmentRelative_AlignmentGeometry_h
+
+#include "FWCore/Utilities/interface/Exception.h"
+
+#include <map>
+#include <set>
+#include <string>
+
+//----------------------------------------------------------------------------------------------------
+
+/**
+ *\brief A structure to hold relevant geometrical information about one detector/sensor.
+ **/
+struct DetGeometry {
+  double z;  ///< z postion at detector centre; mm
+
+  double sx, sy;  ///< detector nominal shift = detector center in global coordinates; in mm
+
+  struct DirectionData {
+    double dx, dy, dz;  ///< x, y and z components of the direction unit vector in global coordinates
+    double s;           ///< projection of (sx, sy) to (dx, dy)
+  };
+
+  std::map<unsigned int, DirectionData> directionData;
+
+  bool isU;  ///< only relevant for strips: true for U detectors, false for V detectors
+             ///< global U, V frame is used - that matches with u, v frame of the 1200 detector
+
+  DetGeometry(double _z = 0., double _sx = 0., double _sy = 0., bool _isU = false)
+      : z(_z), sx(_sx), sy(_sy), isU(_isU) {}
+
+  void setDirection(unsigned int idx, double dx, double dy, double dz) {
+    directionData[idx] = {dx, dy, dz, dx * sx + dy * sy};
+  }
+
+  const DirectionData& getDirectionData(unsigned int idx) const {
+    auto it = directionData.find(idx);
+    if (it == directionData.end())
+      throw cms::Exception("PPS") << "direction index " << idx << " not in the mapping.";
+
+    return it->second;
+  }
+};
+
+//----------------------------------------------------------------------------------------------------
+
+/**
+ * A collection of geometrical information.
+ **/
+class AlignmentGeometry {
+protected:
+  std::map<unsigned int, DetGeometry> sensorGeometry;
+
+public:
+  /// a characteristic z in mm
+  double z0;
+
+  /// puts an element to the map
+  void insert(unsigned int id, const DetGeometry& g);
+
+  /// retrieves sensor geometry
+  const DetGeometry& get(unsigned int id) const;
+
+  const std::map<unsigned int, DetGeometry>& getSensorMap() const { return sensorGeometry; }
+
+  /// returns the number of detectors in the collection
+  unsigned int getNumberOfDetectors() const { return sensorGeometry.size(); }
+
+  /// check whether the sensor Id is valid (present in the map)
+  bool isValidSensorId(unsigned int id) const { return (sensorGeometry.find(id) != sensorGeometry.end()); }
+
+  /// Prints the geometry.
+  void print() const;
+};
+
+#endif

CalibPPS/AlignmentRelative/interface/AlignmentResult.h

@@ -0,0 +1,65 @@
+/****************************************************************************
+* Authors:
+*	Jan Kašpar ([email protected])
+****************************************************************************/
+
+#ifndef CalibPPS_AlignmentRelative_AlignmentResult_h
+#define CalibPPS_AlignmentRelative_AlignmentResult_h
+
+/// \brief Result of CTPPS track-based alignment
+class AlignmentResult {
+protected:
+  /// shift in read-out directions and along beam, mm
+  /// "_unc" denotes the shift uncertainties
+  double sh_r1, sh_r1_unc;
+  double sh_r2, sh_r2_unc;
+  double sh_z, sh_z_unc;
+
+  /// rotation about beam axis, rad
+  double rot_z;
+  double rot_z_unc;
+
+public:
+  AlignmentResult(double _sh_r1 = 0.,
+                  double _sh_r1_e = 0.,
+                  double _sh_r2 = 0.,
+                  double _sh_r2_e = 0.,
+                  double _sh_z = 0.,
+                  double _sh_z_e = 0.,
+                  double _rot_z = 0.,
+                  double _rot_z_e = 0.)
+      : sh_r1(_sh_r1),
+        sh_r1_unc(_sh_r1_e),
+        sh_r2(_sh_r2),
+        sh_r2_unc(_sh_r2_e),
+        sh_z(_sh_z),
+        sh_z_unc(_sh_z_e),
+        rot_z(_rot_z),
+        rot_z_unc(_rot_z_e) {}
+
+  inline double getShR1() const { return sh_r1; }
+  inline void setShR1(const double &v) { sh_r1 = v; }
+
+  inline double getShR1Unc() const { return sh_r1_unc; }
+  inline void setShR1Unc(const double &v) { sh_r1_unc = v; }
+
+  inline double getShR2() const { return sh_r2; }
+  inline void setShR2(const double &v) { sh_r2 = v; }
+
+  inline double getShR2Unc() const { return sh_r2_unc; }
+  inline void setShR2Unc(const double &v) { sh_r2_unc = v; }
+
+  inline double getShZ() const { return sh_z; }
+  inline void setShZ(const double &v) { sh_z = v; }
+
+  inline double getShZUnc() const { return sh_z_unc; }
+  inline void setShZUnc(const double &v) { sh_z_unc = v; }
+
+  inline double getRotZ() const { return rot_z; }
+  inline void setRotZ(const double &v) { rot_z = v; }
+
+  inline double getRotZUnc() const { return rot_z_unc; }
+  inline void setRotZUnc(const double &v) { rot_z_unc = v; }
+};
+
+#endif