diff --git a/CalibPPS/AlignmentRelative/BuildFile.xml b/CalibPPS/AlignmentRelative/BuildFile.xml new file mode 100644 index 0000000000000..de67786cb3d97 --- /dev/null +++ b/CalibPPS/AlignmentRelative/BuildFile.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/CalibPPS/AlignmentRelative/README.md b/CalibPPS/AlignmentRelative/README.md new file mode 100644 index 0000000000000..4f3b9e77c4397 --- /dev/null +++ b/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 + diff --git a/CalibPPS/AlignmentRelative/interface/AlignmentAlgorithm.h b/CalibPPS/AlignmentRelative/interface/AlignmentAlgorithm.h new file mode 100644 index 0000000000000..c5125cc6e23e8 --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/AlignmentAlgorithm.h @@ -0,0 +1,78 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#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 +#include + +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 &, + std::map &results, + TDirectory *dir = nullptr) = 0; + + /// cleans up after processing + virtual void end() = 0; +}; + +#endif diff --git a/CalibPPS/AlignmentRelative/interface/AlignmentConstraint.h b/CalibPPS/AlignmentRelative/interface/AlignmentConstraint.h new file mode 100644 index 0000000000000..f851f83dead28 --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/AlignmentConstraint.h @@ -0,0 +1,29 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#ifndef CalibPPS_AlignmentRelative_AlignmentConstraint_h +#define CalibPPS_AlignmentRelative_AlignmentConstraint_h + +#include + +#include +#include + +/** + *\brief An alignment constraint. + **/ +class AlignmentConstraint { +public: + /// constraint value + double val; + + /// map: AlignmentAlgorithm::QuantityClass -> constraint coefficients + std::map coef; + + /// label of the constraint + std::string name; +}; + +#endif diff --git a/CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h b/CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h new file mode 100644 index 0000000000000..6e454bf1f8d2e --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h @@ -0,0 +1,82 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#ifndef CalibPPS_AlignmentRelative_AlignmentGeometry_h +#define CalibPPS_AlignmentRelative_AlignmentGeometry_h + +#include "FWCore/Utilities/interface/Exception.h" + +#include +#include +#include + +//---------------------------------------------------------------------------------------------------- + +/** + *\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 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 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& 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 diff --git a/CalibPPS/AlignmentRelative/interface/AlignmentResult.h b/CalibPPS/AlignmentRelative/interface/AlignmentResult.h new file mode 100644 index 0000000000000..7d58fc920f2d1 --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/AlignmentResult.h @@ -0,0 +1,65 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#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 diff --git a/CalibPPS/AlignmentRelative/interface/AlignmentTask.h b/CalibPPS/AlignmentRelative/interface/AlignmentTask.h new file mode 100644 index 0000000000000..735b0a15ccfb4 --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/AlignmentTask.h @@ -0,0 +1,133 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#ifndef CalibPPS_AlignmentRelative_AlignmentTask_h +#define CalibPPS_AlignmentRelative_AlignmentTask_h + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h" +class AlignmentConstraint; +class CTPPSGeometry; + +#include + +/** + *\brief Represents an alignment task. + **/ +class AlignmentTask { +public: + // -------------------- config file parameters -------------------- + + /// whether to resolve detector shifts in readout direction(s) + bool resolveShR; + + /// whether to resolve detector shifts in z + bool resolveShZ; + + /// whether to resolve detector rotations around z + bool resolveRotZ; + + /// whether to resolve only 1 rot_z per RP + bool oneRotZPerPot; + + /// whether to apply the constraint mean U = mean V RotZ for strips ("standard" set of constraints only) + bool useEqualMeanUMeanVRotZConstraints; + + /// fixed detectors constraints from config file + edm::ParameterSet fixedDetectorsConstraints; + + /// settings of "standard" constraints from config file + edm::ParameterSet standardConstraints; + + // -------------------- geometry-related members -------------------- + + /// the geometry for this task + AlignmentGeometry geometry; + + /// builds the alignment geometry + static void buildGeometry(const std::vector &rpDecIds, + const std::vector &excludedSensors, + const CTPPSGeometry *, + double z0, + AlignmentGeometry &geometry); + + // -------------------- quantity-class-related members -------------------- + + /// quantity classes + enum QuantityClass { + qcShR1, ///< detector shifts in first readout direction + qcShR2, ///< detector shifts in second readout direction + qcShZ, ///< detector shifts in z + qcRotZ, ///< detector rotations around z + }; + + /// list of quantity classes to be optimized + std::vector quantityClasses; + + /// returns a string tag for the given quantity class + std::string quantityClassTag(QuantityClass) const; + + struct DetIdDirIdxPair { + unsigned int detId; + unsigned int dirIdx; + + bool operator<(const DetIdDirIdxPair &other) const { + if (detId < other.detId) + return true; + if (detId > other.detId) + return false; + if (dirIdx < other.dirIdx) + return true; + + return false; + } + }; + + /// for each quantity class contains mapping (detector id, direction) --> measurement index + std::map> mapMeasurementIndeces; + + /// for each quantity class contains mapping detector id --> quantity index + std::map> mapQuantityIndeces; + + /// builds "mapMatrixIndeces" from "geometry" + void buildIndexMaps(); + + /// returns the number of quantities of the given class + unsigned int measurementsOfClass(QuantityClass) const; + + /// returns the number of quantities of the given class + unsigned int quantitiesOfClass(QuantityClass) const; + + /// returns measurement index (if non-existent, returns -1) + signed int getMeasurementIndex(QuantityClass cl, unsigned int detId, unsigned int dirIdx) const; + + /// returns measurement index (if non-existent, returns -1) + signed int getQuantityIndex(QuantityClass cl, unsigned int detId) const; + + // -------------------- constraint-related members -------------------- + + /// builds a set of fixed-detector constraints + void buildFixedDetectorsConstraints(std::vector &) const; + + /// builds the standard constraints + void buildStandardConstraints(std::vector &) const; + + /// adds constraints such that only 1 rot_z per RP is left + void buildOneRotZPerPotConstraints(std::vector &) const; + + /// adds constraints such that only mean-U and mean-V RotZ are equal for each strip RP + void buildEqualMeanUMeanVRotZConstraints(std::vector &constraints) const; + + // -------------------- constructors -------------------- + + /// dummy constructor (not to be used) + AlignmentTask() {} + + /// normal constructor + AlignmentTask(const edm::ParameterSet &ps); +}; + +#endif diff --git a/CalibPPS/AlignmentRelative/interface/HitCollection.h b/CalibPPS/AlignmentRelative/interface/HitCollection.h new file mode 100644 index 0000000000000..7f303f6e9fc76 --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/HitCollection.h @@ -0,0 +1,39 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#ifndef CalibPPS_AlignmentRelative_HitCollection_h +#define CalibPPS_AlignmentRelative_HitCollection_h + +#include + +//---------------------------------------------------------------------------------------------------- + +struct Hit { + /// sensor id + unsigned int id; + + /// index of read-out direction (valid are: 1 or 2) + unsigned int dirIdx; + + /// measurement position; mm + double position; + + /// measurement position; mm + double sigma; + + /// global z - AlignmentGeometry::z0, mm + double z; + + Hit(unsigned int _id = 0, unsigned int _dirIdx = 0, double _pos = 0, double _sig = 0, double _z = 0) + : id(_id), dirIdx(_dirIdx), position(_pos), sigma(_sig), z(_z) {} +}; + +//---------------------------------------------------------------------------------------------------- + +typedef std::vector HitCollection; + +//---------------------------------------------------------------------------------------------------- + +#endif diff --git a/CalibPPS/AlignmentRelative/interface/IdealResult.h b/CalibPPS/AlignmentRelative/interface/IdealResult.h new file mode 100644 index 0000000000000..629ae187fb014 --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/IdealResult.h @@ -0,0 +1,54 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h" + +#include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionData.h" +#include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsData.h" + +#include "CalibPPS/AlignmentRelative/interface/AlignmentAlgorithm.h" + +#include + +class AlignmentTask; + +/** + *\brief Calculates the ideal result of the StraightTrackAlignment. + **/ +class IdealResult : public AlignmentAlgorithm { +protected: + edm::ESHandle gReal, gMisaligned; + +public: + /// dummy constructor (not to be used) + IdealResult() {} + + /// normal constructor + IdealResult(const edm::ParameterSet &ps, AlignmentTask *_t); + + ~IdealResult() override {} + + std::string getName() override { return "Ideal"; } + + bool hasErrorEstimate() override { return false; } + + void begin(const CTPPSGeometry *geometryReal, const CTPPSGeometry *geometryMisaligned) override; + + void feed(const HitCollection &, const LocalTrackFit &) override {} + + void saveDiagnostics(TDirectory *) override {} + + void analyze() override {} + + unsigned int solve(const std::vector &, + std::map &result, + TDirectory *dir) override; + + void end() override {} +}; diff --git a/CalibPPS/AlignmentRelative/interface/JanAlignmentAlgorithm.h b/CalibPPS/AlignmentRelative/interface/JanAlignmentAlgorithm.h new file mode 100644 index 0000000000000..1d8a678691c1e --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/JanAlignmentAlgorithm.h @@ -0,0 +1,109 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#ifndef CalibPPS_AlignmentRelative_JanAlignmentAlgorithm_h +#define CalibPPS_AlignmentRelative_JanAlignmentAlgorithm_h + +#include "CalibPPS/AlignmentRelative/interface/AlignmentAlgorithm.h" +#include "CalibPPS/AlignmentRelative/interface/SingularMode.h" + +#include "TMatrixD.h" +#include "TVectorD.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TGraph.h" + +#include +#include +#include + +/** + *\brief Jan's alignment algorithm. + **/ +class JanAlignmentAlgorithm : public AlignmentAlgorithm { +public: + /// a scatter plot, with graph and histogram representations + struct ScatterPlot { + TGraph *g; + TH2D *h; + }; + + /// structure holding statistical information for one detector + struct DetStat { + TH1D *m_dist; + TH1D *R_dist; + std::vector coefHist; + std::vector resVsCoef; + std::map, ScatterPlot> resVsCoefRot_perRPSet; + }; + +private: + /// S matrix components + /// indeces correspond to the qClToOpt list + TMatrixD **Sc; + + /// M vector components + /// indeces correspond to the qClToOpt list + TVectorD *Mc; + + /// final S matrix + TMatrixD S; + + /// final M vector + TVectorD M; + + /// eigen values of the S matrix + TVectorD S_eigVal; + + /// matrix of S eigenvectors + TMatrixD S_eigVec; + + /// a list of the singular modes of the S matrix + std::vector singularModes; + + /// whether to stop when singular modes are identified + bool stopOnSingularModes; + + /// normalized eigen value below which the (CS) eigen vectors are considered as weak + double weakLimit; + + /// event count + unsigned int events; + + /// statistical data collection + std::map statistics; + + /// flag whether to build statistical plots + bool buildDiagnosticPlots; + +public: + /// dummy constructor (not to be used) + JanAlignmentAlgorithm() {} + + /// normal constructor + JanAlignmentAlgorithm(const edm::ParameterSet &ps, AlignmentTask *_t); + + ~JanAlignmentAlgorithm() override; + + std::string getName() override { return "Jan"; } + + bool hasErrorEstimate() override { return true; } + + void begin(const CTPPSGeometry *geometryReal, const CTPPSGeometry *geometryMisaligned) override; + + void feed(const HitCollection &, const LocalTrackFit &) override; + + void saveDiagnostics(TDirectory *) override; + + void analyze() override; + + unsigned int solve(const std::vector &, + std::map &results, + TDirectory *dir) override; + + void end() override; +}; + +#endif diff --git a/CalibPPS/AlignmentRelative/interface/LocalTrackFit.h b/CalibPPS/AlignmentRelative/interface/LocalTrackFit.h new file mode 100644 index 0000000000000..6574bd191fecf --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/LocalTrackFit.h @@ -0,0 +1,51 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#ifndef CalibPPS_AlignmentRelative_LocalTrackFit_h +#define CalibPPS_AlignmentRelative_LocalTrackFit_h + +#include "TMath.h" + +/** + *\brief Local (linear) track description (or a fit result). + * Uses global reference system. + **/ +struct LocalTrackFit { + /// the point where intercepts are measured, in mm + double z0; + + /// slopes in rad + double ax, ay; + + /// intercepts in mm + double bx, by; + + /// the number of degrees of freedom + signed int ndf; + + /// the residual sum of squares + double chi_sq; + + LocalTrackFit(double _z0 = 0., + double _ax = 0., + double _ay = 0., + double _bx = 0., + double _by = 0., + unsigned int _ndf = 0, + double _chi_sq = 0.) + : z0(_z0), ax(_ax), ay(_ay), bx(_bx), by(_by), ndf(_ndf), chi_sq(_chi_sq) {} + + double pValue() const { return TMath::Prob(chi_sq, ndf); } + + double chiSqPerNdf() const { return (ndf > 0) ? chi_sq / ndf : 0.; } + + void eval(double z, double &x, double &y) { + double ze = z - z0; + x = ax * ze + bx; + y = ay * ze + by; + } +}; + +#endif diff --git a/CalibPPS/AlignmentRelative/interface/LocalTrackFitter.h b/CalibPPS/AlignmentRelative/interface/LocalTrackFitter.h new file mode 100644 index 0000000000000..dc68f1bd62885 --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/LocalTrackFitter.h @@ -0,0 +1,54 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#ifndef CalibPPS_AlignmentRelative_LocalTrackFitter_h +#define CalibPPS_AlignmentRelative_LocalTrackFitter_h + +#include "CalibPPS/AlignmentRelative/interface/LocalTrackFit.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h" +#include "CalibPPS/AlignmentRelative/interface/HitCollection.h" + +namespace edm { + class ParameterSet; +} + +/** + *\brief Performs straight-line fit and outlier rejection. + **/ +class LocalTrackFitter { +public: + /// dummy constructor (not to be used) + LocalTrackFitter() {} + + /// normal constructor + LocalTrackFitter(const edm::ParameterSet &); + + virtual ~LocalTrackFitter() {} + + /// runs the fit and outlier-removal loop + /// returns true in case of success + bool fit(HitCollection &, const AlignmentGeometry &, LocalTrackFit &) const; + +protected: + /// verbosity level + unsigned int verbosity; + + /// minimum of hits to accept data from a RP + unsigned int minimumHitsPerProjectionPerRP; + + /// hits with higher ratio residual/sigma will be dropped + double maxResidualToSigma; + + /// fits the collection of hits and removes hits with too high residual/sigma ratio + /// \param failed whether the fit has failed + /// \param selectionChanged whether some hits have been removed + void fitAndRemoveOutliers( + HitCollection &, const AlignmentGeometry &, LocalTrackFit &, bool &failed, bool &selectionChanged) const; + + /// removes the hits of pots with too few planes active + void removeInsufficientPots(HitCollection &, bool &selectionChanged) const; +}; + +#endif diff --git a/CalibPPS/AlignmentRelative/interface/SingularMode.h b/CalibPPS/AlignmentRelative/interface/SingularMode.h new file mode 100644 index 0000000000000..84f13e979c8d5 --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/SingularMode.h @@ -0,0 +1,25 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#ifndef CalibPPS_AlignmentRelative_SingularMode_h +#define CalibPPS_AlignmentRelative_SingularMode_h + +#include + +/** + *\brief + **/ +struct SingularMode { + /// eigen value + double val; + + /// eigen vector + TVectorD vec; + + /// index + int idx; +}; + +#endif diff --git a/CalibPPS/AlignmentRelative/interface/StraightTrackAlignment.h b/CalibPPS/AlignmentRelative/interface/StraightTrackAlignment.h new file mode 100644 index 0000000000000..e8b8305a7d9cd --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/StraightTrackAlignment.h @@ -0,0 +1,259 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +* Cristian Baldenegro (crisx.baldenegro@gmail.com) +****************************************************************************/ + +#ifndef CalibPPS_AlignmentRelative_StraightTrackAlignment_h +#define CalibPPS_AlignmentRelative_StraightTrackAlignment_h + +#include +#include +#include + +#include +#include +#include + +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/CTPPSReco/interface/TotemRPRecHit.h" +#include "DataFormats/CTPPSReco/interface/TotemRPUVPattern.h" +#include "DataFormats/CTPPSReco/interface/CTPPSDiamondRecHit.h" +#include "DataFormats/CTPPSReco/interface/CTPPSPixelRecHit.h" +#include "DataFormats/CTPPSReco/interface/CTPPSPixelLocalTrack.h" + +#include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsData.h" + +#include "CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h" +#include "CalibPPS/AlignmentRelative/interface/HitCollection.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentAlgorithm.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentConstraint.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentTask.h" +#include "CalibPPS/AlignmentRelative/interface/LocalTrackFit.h" +#include "CalibPPS/AlignmentRelative/interface/LocalTrackFitter.h" +#include "FWCore/Framework/interface/ESHandle.h" + +namespace edm { + class ParameterSet; + class Event; + class EventSetup; +} // namespace edm + +class TH1D; +class TGraph; + +/** + *\brief Track-based alignment using straight tracks. + **/ +class StraightTrackAlignment { +public: + StraightTrackAlignment(const edm::ParameterSet &); + virtual ~StraightTrackAlignment(); + + virtual void begin(edm::ESHandle hRealAlignment, + edm::ESHandle hRealGeometry, + edm::ESHandle hMisalignedGeometry); + + virtual void processEvent(const edm::EventID &eventId, + const edm::DetSetVector &uvPatternsStrip, + const edm::DetSetVector &hitsDiamond, + const edm::DetSetVector &hitsPixel, + const edm::DetSetVector &tracksPixel); + + /// performs analyses and fill results variable + virtual void finish(); + +protected: + // ---------- input parameters ----------- + + /// verbosity level + unsigned int verbosity; + + /// list of RPs for which the alignment parameters shall be optimized + std::vector rpIds; + + /// list of planes to be excluded from processing + std::vector excludePlanes; + + /// a characteristic z in mm + /// to keep values of z small - this helps the numerical solution + double z0; + + /// the collection of the alignment algorithms + std::vector algorithms; + + /// constraint types + enum ConstraintsType { ctFixedDetectors, ctStandard }; + + /// the chosen type of constraints + ConstraintsType constraintsType; + + /// stops after this event number has been reached + signed int maxEvents; + + // ---------- hit/track selection parameters ---------- + + /// remove events with impossible signatures (i.e. simultaneously top and bottom) + bool removeImpossible; + + /// select only tracks with activity in minimal number of units + unsigned int requireNumberOfUnits; + + /// if a track goes through overlap, select it only if it leaves signal in at least 3 pots + bool requireAtLeast3PotsInOverlap; + + /// if true, only track through vertical-horizontal overlap are seleceted + bool requireOverlap; + + /// whether to cut on chi^2/ndf + bool cutOnChiSqPerNdf; + + /// the value of chi^2/ndf cut threshold + double chiSqPerNdfCut; + + /// cuts on absolute values of the track angle + double maxTrackAx; + double maxTrackAy; + + /// list of RP sets accepted irrespective of the other "require" settings + std::vector > additionalAcceptedRPSets; + + // ---------- output parameters ---------- + + /// file name prefix for result files + std::string fileNamePrefix; + + /// file name prefix for cumulative result files + std::string cumulativeFileNamePrefix; + + /// file name prefix for cumulative expanded result files + std::string expandedFileNamePrefix; + + /// file name prefix for cumulative factored result files + std::string factoredFileNamePrefix; + + /// whether to use long format (many decimal digits) when saving XML files + bool preciseXMLFormat; + + /// whether to save uncertainties in the result XML files + bool saveXMLUncertainties; + + /// whether itermediate results (S, CS matrices) of alignments shall be saved + bool saveIntermediateResults; + + /// the name task data file + std::string taskDataFileName; + + /// the file with task data + TFile *taskDataFile; + + // ---------- internal data members ---------- + + /// the alignment task to be solved + AlignmentTask task; + + /// track fitter + LocalTrackFitter fitter; + + /// (real geometry) alignments before this alignment iteration + CTPPSRPAlignmentCorrectionsData initialAlignments; + + // ---------- diagnostics parameters and plots ---------- + + /// whether to build and save diagnostic plots + bool buildDiagnosticPlots; + + /// file name for some event selection statistics + std::string diagnosticsFile; + + signed int eventsTotal; ///< counter of events + signed int eventsFitted; ///< counter of processed tracks + signed int eventsSelected; ///< counter of processed tracks + + std::map, unsigned long> fittedTracksPerRPSet; ///< counter of fitted tracks in a certain RP set + std::map, unsigned long> + selectedTracksPerRPSet; ///< counter of selected tracks in a certain RP set + + std::map selectedHitsPerPlane; ///< counter of selected hits per plane + + TH1D *fitNdfHist_fitted, *fitNdfHist_selected; ///< fit num. of degrees of freedom histograms for all/selected tracks + TH1D *fitPHist_fitted, *fitPHist_selected; ///< fit p-value histograms for all/selected tracks + TH1D *fitAxHist_fitted, *fitAxHist_selected; ///< fit ax histograms for all/selected tracks + TH1D *fitAyHist_fitted, *fitAyHist_selected; ///< fit ay histograms for all/selected tracks + TH1D *fitBxHist_fitted, *fitBxHist_selected; ///< fit bx histograms for all/selected tracks + TH1D *fitByHist_fitted, *fitByHist_selected; ///< fit by histograms for all/selected tracks + + struct RPSetPlots { + std::string name; + + /// normalised chi^2 histograms for all/selected tracks, in linear/logarithmic scale + TH1D *chisqn_lin_fitted = nullptr, *chisqn_lin_selected = nullptr, *chisqn_log_fitted = nullptr, + *chisqn_log_selected = nullptr; + + /// plots ax vs. ay + TGraph *fitAxVsAyGraph_fitted = nullptr, *fitAxVsAyGraph_selected = nullptr; + + /// plots bx vs. by + TGraph *fitBxVsByGraph_fitted = nullptr, *fitBxVsByGraph_selected = nullptr; + + RPSetPlots() {} + + RPSetPlots(const std::string &_name); + + void free(); + + void write() const; + }; + + /// global (all RP sets) chi^2 histograms + RPSetPlots globalPlots; + + /// chi^2 histograms per RP set + std::map, RPSetPlots> rpSetPlots; + + /// map: detector id --> residua histogram + struct ResiduaHistogramSet { + TH1D *total_fitted, *total_selected; + TGraph *selected_vs_chiSq; + std::map, TH1D *> perRPSet_fitted, perRPSet_selected; + }; + + /// residua histograms + std::map residuaHistograms; + + // ----------- methods ------------ + + /// creates a new residua histogram + TH1D *newResiduaHist(const char *name); + + /// fits the collection of hits and removes hits with too high residual/sigma ratio + /// \param failed whether the fit has failed + /// \param selectionChanged whether some hits have been removed + void fitLocalTrack(HitCollection &, LocalTrackFit &, bool &failed, bool &selectionChanged); + + /// removes the hits of pots with too few planes active + void removeInsufficientPots(HitCollection &, bool &selectionChanged); + + /// builds a selected set of constraints + void buildConstraints(std::vector &); + + /// fills diagnostic (chi^2, residua, ...) histograms + void updateDiagnosticHistograms(const HitCollection &selection, + const std::set &selectedRPs, + const LocalTrackFit &trackFit, + bool trackSelected); + + /// converts a set to string + static std::string setToString(const std::set &); + + /// result pretty printing routines + void printN(const char *str, unsigned int N); + void printLineSeparator(const std::vector > &); + void printQuantitiesLine(const std::vector > &); + void printAlgorithmsLine(const std::vector > &); + + /// saves a ROOT file with diagnostic plots + void saveDiagnostics() const; +}; + +#endif diff --git a/CalibPPS/AlignmentRelative/interface/Utilities.h b/CalibPPS/AlignmentRelative/interface/Utilities.h new file mode 100644 index 0000000000000..73bdda6901fef --- /dev/null +++ b/CalibPPS/AlignmentRelative/interface/Utilities.h @@ -0,0 +1,37 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#ifndef CalibPPS_AlignmentRelative_Utilities_h +#define CalibPPS_AlignmentRelative_Utilities_h + +#include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h" + +#include + +class CTPPSRPAlignmentCorrectionsData; +class AlignmentGeometry; + +extern void printId(unsigned int id); + +extern void print(TMatrixD &m, const char *label = nullptr, bool mathematicaFormat = false); + +/** + * NOTE ON ERROR PROPAGATION + * + * It is not possible to split (and merge again) the experimental errors between the RP and sensor + * contributions. To do so, one would need to keep the entire covariance matrix. Thus, it has been + * decided to save: + * RP errors = the uncertainty of the common shift/rotation + * sensor error = the full experimental uncertainty + * In consequence: RP and sensor errors SHOULD NEVER BE SUMMED! + **/ +extern void factorRPFromSensorCorrections(const CTPPSRPAlignmentCorrectionsData &input, + CTPPSRPAlignmentCorrectionsData &expanded, + CTPPSRPAlignmentCorrectionsData &factored, + const AlignmentGeometry &, + bool equalWeights = false, + unsigned int verbosity = 0); + +#endif diff --git a/CalibPPS/AlignmentRelative/plugins/BuildFile.xml b/CalibPPS/AlignmentRelative/plugins/BuildFile.xml new file mode 100644 index 0000000000000..118547c695728 --- /dev/null +++ b/CalibPPS/AlignmentRelative/plugins/BuildFile.xml @@ -0,0 +1,22 @@ + + + + + + + + + + + + + + + + + + + + + + diff --git a/CalibPPS/AlignmentRelative/plugins/PPSFastLocalSimulation.cc b/CalibPPS/AlignmentRelative/plugins/PPSFastLocalSimulation.cc new file mode 100644 index 0000000000000..7dd378ac75c1e --- /dev/null +++ b/CalibPPS/AlignmentRelative/plugins/PPSFastLocalSimulation.cc @@ -0,0 +1,440 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/ESGetToken.h" +#include "FWCore/Utilities/interface/RandomNumberGenerator.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "CLHEP/Random/RandomEngine.h" +#include "CLHEP/Random/RandGauss.h" +#include "CLHEP/Random/RandExponential.h" + +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" + +#include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h" +#include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" + +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/CTPPSReco/interface/TotemRPRecHit.h" +#include "DataFormats/CTPPSReco/interface/CTPPSDiamondRecHit.h" +#include "DataFormats/CTPPSReco/interface/CTPPSPixelRecHit.h" + +#include "Geometry/VeryForwardRPTopology/interface/RPTopology.h" +#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h" +#include "Geometry/Records/interface/VeryForwardMisalignedGeometryRecord.h" + +#include "CalibPPS/AlignmentRelative/interface/Utilities.h" + +#include "TMath.h" +#include "TMatrixD.h" +#include "TVectorD.h" + +/** + *\brief Fast (no G4) proton simulation in within one station. + * Uses misaligned geometry. + */ +class PPSFastLocalSimulation : public edm::stream::EDProducer<> { +public: + PPSFastLocalSimulation(const edm::ParameterSet &); + ~PPSFastLocalSimulation() override; + +protected: + /// verbosity level + unsigned int verbosity_; + + /// whether a HepMC description of the proton shall be saved in the event + bool makeHepMC_; + + /// whether the hits of the proton shall be calculated and saved + bool makeHits_; + + /// the list of RPs to simulate + std::vector RPs_; + + /// number of particles to generate per event + unsigned int particlesPerEvent_; + + /// particle energy and momentum + double particle_E_, particle_p_; + + /// the "origin" of tracks, in mm + double z0_; + + /// whether measurement values shall be rounded to the nearest strip + bool roundToPitch_; + + /// in mm + double pitchStrips_, pitchDiamonds_, pitchPixels_; + + /// size of insensitive margin at sensor's edge facing the beam, in mm + double insensitiveMarginStrips_; + + struct Distribution { + enum Type { dtBox, dtGauss, dtGaussLimit } type_; + double x_mean_, x_width_, x_min_, x_max_; + double y_mean_, y_width_, y_min_, y_max_; + + Distribution(const edm::ParameterSet &); + + void Generate(CLHEP::HepRandomEngine &rndEng, double &x, double &y); + }; + + /// position parameters in mm + Distribution position_dist_; + + /// angular parameters in rad + Distribution angular_dist_; + + //---------- internal parameters ---------- + + /// v position of strip 0, in mm + double stripZeroPosition_; + + edm::ESGetToken esTokenGeometry_; + + void GenerateTrack(unsigned int pi, + CLHEP::HepRandomEngine &rndEng, + HepMC::GenEvent *gEv, + std::unique_ptr> &stripHitColl, + std::unique_ptr> &diamondHitColl, + std::unique_ptr> &pixelHitColl, + const CTPPSGeometry &geometry); + + //---------- framework methods ---------- + + void produce(edm::Event &, const edm::EventSetup &) override; +}; + +//---------------------------------------------------------------------------------------------------- + +using namespace edm; +using namespace std; +using namespace CLHEP; +using namespace HepMC; + +//---------------------------------------------------------------------------------------------------- + +PPSFastLocalSimulation::Distribution::Distribution(const edm::ParameterSet &ps) { + // get type + string typeName = ps.getParameter("type"); + if (!typeName.compare("box")) + type_ = dtBox; + else if (!typeName.compare("gauss")) + type_ = dtGauss; + else if (!typeName.compare("gauss-limit")) + type_ = dtGaussLimit; + else + throw cms::Exception("PPS") << "Unknown distribution type `" << typeName << "'."; + + x_mean_ = ps.getParameter("x_mean"); + x_width_ = ps.getParameter("x_width"); + x_min_ = ps.getParameter("x_min"); + x_max_ = ps.getParameter("x_max"); + + y_mean_ = ps.getParameter("y_mean"); + y_width_ = ps.getParameter("y_width"); + y_min_ = ps.getParameter("y_min"); + y_max_ = ps.getParameter("y_max"); +} + +//---------------------------------------------------------------------------------------------------- + +void PPSFastLocalSimulation::Distribution::Generate(CLHEP::HepRandomEngine &rndEng, double &x, double &y) { + switch (type_) { + case dtBox: + x = x_mean_ + x_width_ * (rndEng.flat() - 0.5); + y = y_mean_ + y_width_ * (rndEng.flat() - 0.5); + break; + + case dtGauss: + x = x_mean_ + RandGauss::shoot(&rndEng) * x_width_; + y = y_mean_ + RandGauss::shoot(&rndEng) * y_width_; + break; + + case dtGaussLimit: { + const double u_x = rndEng.flat(), u_y = rndEng.flat(); + + const double cdf_x_min = (1. + TMath::Erf((x_min_ - x_mean_) / x_width_ / sqrt(2.))) / 2.; + const double cdf_x_max = (1. + TMath::Erf((x_max_ - x_mean_) / x_width_ / sqrt(2.))) / 2.; + const double a_x = cdf_x_max - cdf_x_min, b_x = cdf_x_min; + + const double cdf_y_min = (1. + TMath::Erf((y_min_ - y_mean_) / y_width_ / sqrt(2.))) / 2.; + const double cdf_y_max = (1. + TMath::Erf((y_max_ - y_mean_) / y_width_ / sqrt(2.))) / 2.; + const double a_y = cdf_y_max - cdf_y_min, b_y = cdf_y_min; + + x = x_mean_ + x_width_ * sqrt(2.) * TMath::ErfInverse(2. * (a_x * u_x + b_x) - 1.); + y = y_mean_ + y_width_ * sqrt(2.) * TMath::ErfInverse(2. * (a_y * u_y + b_y) - 1.); + } + + break; + + default: + x = y = 0.; + } +} + +//---------------------------------------------------------------------------------------------------- + +PPSFastLocalSimulation::PPSFastLocalSimulation(const edm::ParameterSet &ps) + : verbosity_(ps.getUntrackedParameter("verbosity", 0)), + + makeHepMC_(ps.getParameter("makeHepMC")), + makeHits_(ps.getParameter("makeHits")), + + RPs_(ps.getParameter>("RPs")), + + particlesPerEvent_(ps.getParameter("particlesPerEvent")), + particle_E_(ps.getParameter("particle_E")), + particle_p_(ps.getParameter("particle_p")), + z0_(ps.getParameter("z0")), + + roundToPitch_(ps.getParameter("roundToPitch")), + pitchStrips_(ps.getParameter("pitchStrips")), + pitchDiamonds_(ps.getParameter("pitchDiamonds")), + pitchPixels_(ps.getParameter("pitchPixels")), + + insensitiveMarginStrips_(ps.getParameter("insensitiveMarginStrips")), + + position_dist_(ps.getParameterSet("position_distribution")), + angular_dist_(ps.getParameterSet("angular_distribution")), + + esTokenGeometry_(esConsumes()) { + // v position of strip 0 + stripZeroPosition_ = RPTopology::last_strip_to_border_dist_ + (RPTopology::no_of_strips_ - 1) * RPTopology::pitch_ - + RPTopology::y_width_ / 2.; + + // register the output + if (makeHepMC_) + produces(); + + if (makeHits_) { + produces>(); + produces>(); + produces>(); + } +} + +//---------------------------------------------------------------------------------------------------- + +PPSFastLocalSimulation::~PPSFastLocalSimulation() {} + +//---------------------------------------------------------------------------------------------------- + +void PPSFastLocalSimulation::GenerateTrack(unsigned int idx, + CLHEP::HepRandomEngine &rndEng, + HepMC::GenEvent *gEv, + unique_ptr> &stripHitColl, + unique_ptr> &diamondHitColl, + unique_ptr> &pixelHitColl, + const CTPPSGeometry &geometry) { + // generate track + double bx = 0., by = 0., ax = 0., ay = 0.; + position_dist_.Generate(rndEng, bx, by); + angular_dist_.Generate(rndEng, ax, ay); + + if (verbosity_ > 5) + printf("\tax = %.3f mrad, bx = %.3f mm, ay = %.3f mrad, by = %.3f mm, z0 = %.3f m\n", + ax * 1E3, + bx, + ay * 1E3, + by, + z0_ * 1E-3); + + // add HepMC track description + if (makeHepMC_) { + GenVertex *gVx = new GenVertex(HepMC::FourVector(bx, by, z0_, 0.)); + gEv->add_vertex(gVx); + + GenParticle *gPe; + double az = sqrt(1. - ax * ax - ay * ay); + gPe = new GenParticle(HepMC::FourVector(particle_p_ * ax, particle_p_ * ay, particle_p_ * az, particle_E_), + 2212, + 1); // add a proton in final state + gPe->suggest_barcode(idx + 1); + gVx->add_particle_out(gPe); + } + + if (makeHits_) { + // check all sensors known to geometry + for (CTPPSGeometry::mapType::const_iterator it = geometry.beginSensor(); it != geometry.endSensor(); ++it) { + // get RP decimal id + CTPPSDetId detId(it->first); + unsigned int decRPId = detId.arm() * 100 + detId.station() * 10 + detId.rp(); + + // stop if the RP is not selected + if (find(RPs_.begin(), RPs_.end(), decRPId) == RPs_.end()) + continue; + + // keep only 1 diamond channel to represent 1 plane + if (detId.subdetId() == CTPPSDetId::sdTimingDiamond) { + CTPPSDiamondDetId channelId(it->first); + if (channelId.channel() != 0) + continue; + } + + if (verbosity_ > 5) { + printf(" "); + printId(it->first); + printf(": "); + } + + // determine the track impact point (in global coordinates) + // !! this assumes that local axes (1, 0, 0) and (0, 1, 0) describe the sensor surface + const auto gl_o = geometry.localToGlobal(detId, CTPPSGeometry::Vector(0, 0, 0)); + const auto gl_a1 = geometry.localToGlobal(detId, CTPPSGeometry::Vector(1, 0, 0)) - gl_o; + const auto gl_a2 = geometry.localToGlobal(detId, CTPPSGeometry::Vector(0, 1, 0)) - gl_o; + + TMatrixD A(3, 3); + TVectorD B(3); + A(0, 0) = ax; + A(0, 1) = -gl_a1.x(); + A(0, 2) = -gl_a2.x(); + B(0) = gl_o.x() - bx; + A(1, 0) = ay; + A(1, 1) = -gl_a1.y(); + A(1, 2) = -gl_a2.y(); + B(1) = gl_o.y() - by; + A(2, 0) = 1.; + A(2, 1) = -gl_a1.z(); + A(2, 2) = -gl_a2.z(); + B(2) = gl_o.z() - z0_; + TMatrixD Ai(3, 3); + Ai = A.Invert(); + TVectorD P(3); + P = Ai * B; + + double de_z = P(0); + CTPPSGeometry::Vector h_glo(ax * de_z + bx, ay * de_z + by, de_z + z0_); + + // hit in local coordinates + CTPPSGeometry::Vector h_loc = geometry.globalToLocal(detId, h_glo); + + // strips + if (detId.subdetId() == CTPPSDetId::sdTrackingStrip) { + double u = h_loc.x(); + double v = h_loc.y(); + + if (verbosity_ > 5) + printf(" u=%+8.4f, v=%+8.4f", u, v); + + // is it within detector? + if (!RPTopology::IsHit(u, v, insensitiveMarginStrips_)) { + if (verbosity_ > 5) + printf(" | no hit\n"); + continue; + } + + // round the measurement + if (roundToPitch_) { + double m = stripZeroPosition_ - v; + signed int strip = (int)floor(m / pitchStrips_ + 0.5); + + v = stripZeroPosition_ - pitchStrips_ * strip; + + if (verbosity_ > 5) + printf(" | strip=%+4i", strip); + } + + double sigma = pitchStrips_ / sqrt(12.); + + if (verbosity_ > 5) + printf(" | m=%+8.4f, sigma=%+8.4f\n", v, sigma); + + DetSet &hits = stripHitColl->find_or_insert(detId); + hits.push_back(TotemRPRecHit(v, sigma)); + } + + // diamonds + if (detId.subdetId() == CTPPSDetId::sdTimingDiamond) { + if (roundToPitch_) { + h_loc.SetX(pitchDiamonds_ * floor(h_loc.x() / pitchDiamonds_ + 0.5)); + } + + if (verbosity_ > 5) + printf(" m = %.3f\n", h_loc.x()); + + const double width = pitchDiamonds_; + + DetSet &hits = diamondHitColl->find_or_insert(detId); + HPTDCErrorFlags flags; + hits.push_back(CTPPSDiamondRecHit(h_loc.x(), width, 0., 0., 0., 0., 0., 0., 0., 0, flags, false)); + } + + // pixels + if (detId.subdetId() == CTPPSDetId::sdTrackingPixel) { + if (roundToPitch_) { + h_loc.SetX(pitchPixels_ * floor(h_loc.x() / pitchPixels_ + 0.5)); + h_loc.SetY(pitchPixels_ * floor(h_loc.y() / pitchPixels_ + 0.5)); + } + + if (verbosity_ > 5) + printf(" m1 = %.3f, m2 = %.3f\n", h_loc.x(), h_loc.y()); + + const double sigma = pitchPixels_ / sqrt(12.); + + const LocalPoint lp(h_loc.x(), h_loc.y(), h_loc.z()); + const LocalError le(sigma, 0., sigma); + + DetSet &hits = pixelHitColl->find_or_insert(detId); + hits.push_back(CTPPSPixelRecHit(lp, le)); + } + } + } +} + +//---------------------------------------------------------------------------------------------------- + +void PPSFastLocalSimulation::produce(edm::Event &event, const edm::EventSetup &es) { + if (verbosity_ > 2) + printf(">> PPSFastLocalSimulation::produce > event %llu\n", event.id().event()); + + Service rng; + HepRandomEngine &rndEng = rng->getEngine(event.streamID()); + + if (verbosity_ > 5) + printf("\tseed = %li\n", rndEng.getSeed()); + + // get geometry + auto const &geometry = es.getData(esTokenGeometry_); + + // initialize products + GenEvent *gEv = new GenEvent(); + gEv->set_event_number(event.id().event()); + + unique_ptr> stripHitColl(new DetSetVector()); + unique_ptr> diamondHitColl(new DetSetVector()); + unique_ptr> pixelHitColl(new DetSetVector()); + + // run particle loop + for (unsigned int pi = 0; pi < particlesPerEvent_; pi++) { + if (verbosity_ > 5) + printf(" generating track %u\n", pi); + + GenerateTrack(pi, rndEng, gEv, stripHitColl, diamondHitColl, pixelHitColl, geometry); + } + + // save products + if (makeHepMC_) { + unique_ptr hepMCoutput(new HepMCProduct()); + hepMCoutput->addHepMCData(gEv); + event.put(move(hepMCoutput)); + } + + if (makeHits_) { + event.put(move(stripHitColl)); + event.put(move(diamondHitColl)); + event.put(move(pixelHitColl)); + } +} + +//---------------------------------------------------------------------------------------------------- + +DEFINE_FWK_MODULE(PPSFastLocalSimulation); diff --git a/CalibPPS/AlignmentRelative/plugins/PPSModifySingularModes.cc b/CalibPPS/AlignmentRelative/plugins/PPSModifySingularModes.cc new file mode 100644 index 0000000000000..b2f2b1dc5c59c --- /dev/null +++ b/CalibPPS/AlignmentRelative/plugins/PPSModifySingularModes.cc @@ -0,0 +1,140 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "FWCore/Framework/interface/stream/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/ParameterSet/interface/FileInPath.h" + +#include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionData.h" +#include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsData.h" + +#include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h" +#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h" +#include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsMethods.h" + +#include "CalibPPS/AlignmentRelative/interface/AlignmentTask.h" +#include "CalibPPS/AlignmentRelative/interface/Utilities.h" + +/** + *\brief Modifies the alignment modes unconstrained by the track-based alignment. + **/ +class PPSModifySingularModes : public edm::stream::EDAnalyzer<> { +public: + PPSModifySingularModes(const edm::ParameterSet &ps); + +private: + edm::ParameterSet ps_; + + edm::ESGetToken tokenRealGeometry_; + + void beginRun(edm::Run const &, edm::EventSetup const &) override; + + void analyze(const edm::Event &e, const edm::EventSetup &es) override {} +}; + +using namespace std; +using namespace edm; + +//---------------------------------------------------------------------------------------------------- + +PPSModifySingularModes::PPSModifySingularModes(const ParameterSet &ps) + : ps_(ps), tokenRealGeometry_(esConsumes()) {} + +//---------------------------------------------------------------------------------------------------- + +void PPSModifySingularModes::beginRun(edm::Run const &, edm::EventSetup const &es) { + // get config parameters + const double z1 = ps_.getUntrackedParameter("z1"); + const double z2 = ps_.getUntrackedParameter("z2"); + const double de_x1 = ps_.getUntrackedParameter("de_x1"); + const double de_x2 = ps_.getUntrackedParameter("de_x2"); + const double de_y1 = ps_.getUntrackedParameter("de_y1"); + const double de_y2 = ps_.getUntrackedParameter("de_y2"); + const double de_rho1 = ps_.getUntrackedParameter("de_rho1"); + const double de_rho2 = ps_.getUntrackedParameter("de_rho2"); + + FileInPath inputFileInPath(ps_.getUntrackedParameter("inputFile")); + const string inputFile = inputFileInPath.fullPath(); + const string outputFile = ps_.getUntrackedParameter("outputFile"); + + // validate config parameters + if (z1 == z2) + throw cms::Exception("PPS") << "z1 equals z2"; + + // calculate slopes and intercepts + const double a_x = (de_x2 - de_x1) / (z2 - z1), b_x = de_x1 - a_x * z1; + const double a_y = (de_y2 - de_y1) / (z2 - z1), b_y = de_y1 - a_y * z1; + const double a_rho = (de_rho2 - de_rho1) / (z2 - z1), b_rho = de_rho1 - a_rho * z1; + + // get geometry + const auto &geometry = es.getData(tokenRealGeometry_); + + // get input alignments + CTPPSRPAlignmentCorrectionsDataSequence inputSequence = CTPPSRPAlignmentCorrectionsMethods::loadFromXML(inputFile); + const auto &input = inputSequence.begin()->second; + + // modify the singular modes + CTPPSRPAlignmentCorrectionsData output = input; + + for (auto &it : input.getSensorMap()) { + const auto &sensorId = it.first; + + const auto &c = geometry.sensorTranslation(sensorId); + + // pixels cannot be described by one single value of z, but no better approxiamtion can be made + const double z = c.z(); + + double de_ShX = a_x * z + b_x; + double de_ShY = a_y * z + b_y; + const double de_RotZ = a_rho * z + b_rho; + + // add the effect of global rotation (about origin, not sensor centre) + de_ShX -= +de_RotZ * (c.y() + de_ShY); + de_ShY -= -de_RotZ * (c.x() + de_ShX); + + CTPPSRPAlignmentCorrectionData d = it.second; + d.setShX(d.getShX() + de_ShX); + d.setShY(d.getShY() + de_ShY); + d.setRotZ(d.getRotZ() + de_RotZ); + + output.setSensorCorrection(sensorId, d); + } + + // build list of RPs + vector rps; + unsigned int last_rp = 123456; + for (auto &it : input.getSensorMap()) { + CTPPSDetId senId(it.first); + unsigned int rpDecId = senId.arm() * 100 + senId.station() * 10 + senId.rp(); + + if (last_rp != rpDecId) { + rps.push_back(rpDecId); + last_rp = rpDecId; + } + } + + // build alignment geometry (needed for the factorisation below) + AlignmentGeometry alignmentGeometry; + vector excludePlanes; + AlignmentTask::buildGeometry(rps, excludePlanes, &geometry, 0., alignmentGeometry); + + // factorise output + CTPPSRPAlignmentCorrectionsData outputExpanded; + CTPPSRPAlignmentCorrectionsData outputFactored; + + const bool equalWeights = false; + factorRPFromSensorCorrections(output, outputExpanded, outputFactored, alignmentGeometry, equalWeights, 1); + + // save output + CTPPSRPAlignmentCorrectionsMethods::writeToXML(outputFactored, outputFile, false, false, true, true, true, true); +} + +//---------------------------------------------------------------------------------------------------- + +DEFINE_FWK_MODULE(PPSModifySingularModes); diff --git a/CalibPPS/AlignmentRelative/plugins/PPSStraightTrackAligner.cc b/CalibPPS/AlignmentRelative/plugins/PPSStraightTrackAligner.cc new file mode 100644 index 0000000000000..8aebb7a350763 --- /dev/null +++ b/CalibPPS/AlignmentRelative/plugins/PPSStraightTrackAligner.cc @@ -0,0 +1,160 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +* Cristian Baldenegro (crisx.baldenegro@gmail.com) +****************************************************************************/ + +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/ESWatcher.h" +#include "FWCore/Utilities/interface/ESGetToken.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include "CalibPPS/AlignmentRelative/interface/StraightTrackAlignment.h" + +#include "CondFormats/AlignmentRecord/interface/RPRealAlignmentRecord.h" +#include "Geometry/Records/interface/VeryForwardMisalignedGeometryRecord.h" +#include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h" +#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h" + +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/CTPPSReco/interface/TotemRPRecHit.h" +#include "DataFormats/CTPPSReco/interface/TotemRPUVPattern.h" +#include "DataFormats/CTPPSReco/interface/CTPPSDiamondRecHit.h" +#include "DataFormats/CTPPSReco/interface/CTPPSPixelRecHit.h" +#include "DataFormats/CTPPSReco/interface/CTPPSPixelLocalTrack.h" + +/** + *\brief An EDAnalyzer that runs StraightTrackAlignment. + **/ +class PPSStraightTrackAligner : public edm::one::EDAnalyzer<> { +public: + PPSStraightTrackAligner(const edm::ParameterSet &ps); + ~PPSStraightTrackAligner() override {} + +private: + unsigned int verbosity_; + + edm::InputTag tagUVPatternsStrip_; + edm::InputTag tagDiamondHits_; + edm::InputTag tagPixelHits_; + edm::InputTag tagPixelLocalTracks_; + + edm::EDGetTokenT> tokenUVPatternsStrip_; + edm::EDGetTokenT> tokenDiamondHits_; + edm::EDGetTokenT> tokenPixelHits_; + edm::EDGetTokenT> tokenPixelLocalTracks_; + + edm::ESWatcher geometryWatcher_; + + edm::ESGetToken esTokenRealGeometry_; + edm::ESGetToken esTokenMisalignedGeometry_; + edm::ESGetToken esTokenRealAlignment_; + + bool worker_initialized_; + StraightTrackAlignment worker_; + + void analyze(const edm::Event &e, const edm::EventSetup &es) override; + + void endJob() override; +}; + +using namespace std; +using namespace edm; + +//---------------------------------------------------------------------------------------------------- + +PPSStraightTrackAligner::PPSStraightTrackAligner(const ParameterSet &ps) + : verbosity_(ps.getUntrackedParameter("verbosity", 0)), + + tagUVPatternsStrip_(ps.getParameter("tagUVPatternsStrip")), + tagDiamondHits_(ps.getParameter("tagDiamondHits")), + tagPixelHits_(ps.getParameter("tagPixelHits")), + tagPixelLocalTracks_(ps.getParameter("tagPixelLocalTracks")), + + tokenUVPatternsStrip_(consumes>(tagUVPatternsStrip_)), + tokenDiamondHits_(consumes>(tagDiamondHits_)), + tokenPixelHits_(consumes>(tagPixelHits_)), + tokenPixelLocalTracks_(consumes>(tagPixelLocalTracks_)), + + esTokenRealGeometry_(esConsumes()), + esTokenMisalignedGeometry_(esConsumes()), + esTokenRealAlignment_(esConsumes()), + + worker_initialized_(false), + worker_(ps) { + if (!tagPixelHits_.label().empty() && !tagPixelLocalTracks_.label().empty()) + LogWarning("PPS") + << "Both tagPixelHits and tagPixelLocalTracks are not empty - most likely this not what you want."; +} + +//---------------------------------------------------------------------------------------------------- + +void PPSStraightTrackAligner::analyze(const edm::Event &event, const edm::EventSetup &es) { + // check if geometry hasn't changed + if (geometryWatcher_.check(es)) { + if (worker_initialized_) + throw cms::Exception("PPS") << "PPSStraightTrackAligner can't cope with changing geometry - change in event " + << event.id() << endl; + } + + // check if worker already initialised + if (!worker_initialized_) { + auto hRealGeometry = es.getHandle(esTokenRealGeometry_); + auto hMisalignedGeometry = es.getHandle(esTokenMisalignedGeometry_); + auto hRealAlignment = es.getHandle(esTokenRealAlignment_); + + worker_.begin(hRealAlignment, hRealGeometry, hMisalignedGeometry); + worker_initialized_ = true; + } + + // get input + DetSetVector defaultStripUVPatterns; + const DetSetVector *pStripUVPatterns = &defaultStripUVPatterns; + Handle> inputStripUVPatterns; + if (!tagUVPatternsStrip_.label().empty()) { + event.getByToken(tokenUVPatternsStrip_, inputStripUVPatterns); + pStripUVPatterns = &(*inputStripUVPatterns); + } + + DetSetVector defaultDiamondHits; + const DetSetVector *pDiamondHits = &defaultDiamondHits; + Handle> inputDiamondHits; + if (!tagDiamondHits_.label().empty()) { + event.getByToken(tokenDiamondHits_, inputDiamondHits); + pDiamondHits = &(*inputDiamondHits); + } + + DetSetVector defaultPixelHits; + const DetSetVector *pPixelHits = &defaultPixelHits; + Handle> inputPixelHits; + if (!tagPixelHits_.label().empty()) { + event.getByToken(tokenPixelHits_, inputPixelHits); + pPixelHits = &(*inputPixelHits); + } + + DetSetVector defaultPixelLocalTracks; + const DetSetVector *pPixelLocalTracks = &defaultPixelLocalTracks; + Handle> inputPixelLocalTracks; + if (!tagPixelLocalTracks_.label().empty()) { + event.getByToken(tokenPixelLocalTracks_, inputPixelLocalTracks); + pPixelLocalTracks = &(*inputPixelLocalTracks); + } + + // feed worker + worker_.processEvent(event.id(), *pStripUVPatterns, *pDiamondHits, *pPixelHits, *pPixelLocalTracks); +} + +//---------------------------------------------------------------------------------------------------- + +void PPSStraightTrackAligner::endJob() { + if (worker_initialized_) + worker_.finish(); + else + throw cms::Exception("PPS") << "worker not initialized." << endl; +} + +DEFINE_FWK_MODULE(PPSStraightTrackAligner); diff --git a/CalibPPS/AlignmentRelative/python/ppsFastLocalSimulation_cfi.py b/CalibPPS/AlignmentRelative/python/ppsFastLocalSimulation_cfi.py new file mode 100644 index 0000000000000..e46a7d7c35978 --- /dev/null +++ b/CalibPPS/AlignmentRelative/python/ppsFastLocalSimulation_cfi.py @@ -0,0 +1,52 @@ +import FWCore.ParameterSet.Config as cms + +ppsFastLocalSimulation = cms.EDProducer("PPSFastLocalSimulation", + verbosity = cms.untracked.uint32(0), + + makeHepMC = cms.bool(False), + makeHits = cms.bool(True), + + particlesPerEvent = cms.uint32(1), + + particle_p = cms.double(6500), # in GeV + particle_E = cms.double(6500), # in GeV + + z0 = cms.double(214500), + + position_distribution = cms.PSet( + type = cms.string("box"), + x_mean = cms.double(5.0), #in mm + x_width = cms.double(10.0), + x_min = cms.double(0.0), + x_max = cms.double(0.0), + + y_mean = cms.double(0.0), + y_width = cms.double(20.0), + y_min = cms.double(0.0), + y_max = cms.double(0.0) + ), + + angular_distribution = cms.PSet( + type = cms.string("gauss"), + x_mean = cms.double(0.0), #in rad + x_width = cms.double(100E-6), + x_min = cms.double(0E-6), + x_max = cms.double(0E-6), + + y_mean = cms.double(0.0), + y_width = cms.double(100E-6), + y_min = cms.double(0E-6), + y_max = cms.double(0E-6) + ), + + #RPs = cms.vuint32(120, 121, 122, 123, 124, 125), + RPs = cms.vuint32(103, 116, 123), + + roundToPitch = cms.bool(True), + + pitchStrips = cms.double(66E-3), # mm + pitchDiamonds = cms.double(200E-3), # mm + pitchPixels = cms.double(30E-3), # mm + + insensitiveMarginStrips = cms.double(34E-3), #mm, save value as RPActiveEdgePosition in SimTotem/RPDigiProducer/python/RPSiDetConf_cfi.py +) diff --git a/CalibPPS/AlignmentRelative/python/ppsModifySingularModes_cfi.py b/CalibPPS/AlignmentRelative/python/ppsModifySingularModes_cfi.py new file mode 100644 index 0000000000000..b86e061126bb9 --- /dev/null +++ b/CalibPPS/AlignmentRelative/python/ppsModifySingularModes_cfi.py @@ -0,0 +1,20 @@ +import FWCore.ParameterSet.Config as cms + +ppsModifySingularModes = cms.EDAnalyzer("PPSModifySingularModes", + inputFile = cms.untracked.string(""), + outputFile = cms.untracked.string(""), + + # mm + z1 = cms.untracked.double(0), + z2 = cms.untracked.double(0), + + # mm + de_x1 = cms.untracked.double(0), + de_y1 = cms.untracked.double(0), + de_x2 = cms.untracked.double(0), + de_y2 = cms.untracked.double(0), + + # rad + de_rho1 = cms.untracked.double(0), + de_rho2 = cms.untracked.double(0), +) diff --git a/CalibPPS/AlignmentRelative/python/ppsStraightTrackAligner_cfi.py b/CalibPPS/AlignmentRelative/python/ppsStraightTrackAligner_cfi.py new file mode 100644 index 0000000000000..4aad7ae75b4e5 --- /dev/null +++ b/CalibPPS/AlignmentRelative/python/ppsStraightTrackAligner_cfi.py @@ -0,0 +1,131 @@ +import FWCore.ParameterSet.Config as cms + +ppsStraightTrackAligner = cms.EDAnalyzer("PPSStraightTrackAligner", + verbosity = cms.untracked.uint32(1), + + # ---------- input and event selection ---------- + + tagUVPatternsStrip = cms.InputTag("totemRPUVPatternFinder"), + tagDiamondHits = cms.InputTag("ctppsFastLocalSimulation"), + tagPixelHits = cms.InputTag("ctppsFastLocalSimulation"), + tagPixelLocalTracks = cms.InputTag(""), + + # list of RPs for which the alignment parameters shall be optimized + rpIds = cms.vuint32(), + + # list of planes to be excluded from processing + excludePlanes = cms.vuint32(), + + # maximum number of selected events + maxEvents = cms.int32(-1), # -1 means unlimited + + + # ---------- event selection ---------- + + # parameters of hit outlier removal + maxResidualToSigma = cms.double(3), + minimumHitsPerProjectionPerRP = cms.uint32(4), + + # skip events with hits in both top and bottom RPs + removeImpossible = cms.bool(True), + + # minimum required number of units active + requireNumberOfUnits = cms.uint32(2), + + # require combination of top+horizontal or bottom+horizontal RPs + requireOverlap = cms.bool(False), + + # require at least 3 RPs active when track in the horizontal-vertical overlap + requireAtLeast3PotsInOverlap = cms.bool(True), + + # list of RP sets that are accepted irrespective of the "require" settings + # the sets should be separated by ";" + # within each set, RP ids are separated by "," + # example: "103,104;120,121" + additionalAcceptedRPSets = cms.string(""), + + # track fit quality requirements + cutOnChiSqPerNdf = cms.bool(True), + chiSqPerNdfCut = cms.double(10), + + # track angular requirements + maxTrackAx = cms.double(1E6), + maxTrackAy = cms.double(1E6), + + + # ---------- constraints ---------- + + # choices: fixedDetectors, standard + constraintsType = cms.string("standard"), + + oneRotZPerPot = cms.bool(False), + useEqualMeanUMeanVRotZConstraints = cms.bool(True), + + # values in um and m rad + fixedDetectorsConstraints = cms.PSet( + ShR = cms.PSet( + ids = cms.vuint32(1220, 1221, 1228, 1229), + values = cms.vdouble(0, 0, 0, 0), + ), + ShZ = cms.PSet( + ids = cms.vuint32(1200, 1201, 1208, 1209), + values = cms.vdouble(0, 0, 0, 0), + ), + RotZ = cms.PSet( + ids = cms.vuint32(1200, 1201), + values = cms.vdouble(0, 0), + ), + RPShZ = cms.PSet( + ids = cms.vuint32(1200), # number of any plane in the chosen RP + values = cms.vdouble(0), + ), + ), + + standardConstraints = cms.PSet( + units = cms.vuint32(1, 21) + ), + + + # ---------- solution parameters ---------- + + # a characteristic z in mm + z0 = cms.double(0.0), + + # what to be resolved + resolveShR = cms.bool(True), + resolveShZ = cms.bool(False), + resolveRotZ = cms.bool(False), + + # suitable value for station alignment + singularLimit = cms.double(1E-8), + + + # ---------- algorithm configuration ---------- + + # available algorithms: Ideal, Jan + algorithms = cms.vstring("Jan"), + + IdealResult = cms.PSet( + ), + + JanAlignmentAlgorithm = cms.PSet( + weakLimit = cms.double(1E-6), + stopOnSingularModes = cms.bool(True), + buildDiagnosticPlots = cms.bool(True), + ), + + + # ---------- output configuration ---------- + + saveIntermediateResults = cms.bool(True), + taskDataFileName = cms.string(''), + + buildDiagnosticPlots = cms.bool(True), + diagnosticsFile = cms.string(''), + + fileNamePrefix = cms.string('results_iteration_'), + expandedFileNamePrefix = cms.string('results_cumulative_expanded_'), + factoredFileNamePrefix = cms.string('results_cumulative_factored_'), + preciseXMLFormat = cms.bool(False), + saveXMLUncertainties = cms.bool(False) +) diff --git a/CalibPPS/AlignmentRelative/src/AlignmentAlgorithm.cc b/CalibPPS/AlignmentRelative/src/AlignmentAlgorithm.cc new file mode 100644 index 0000000000000..0a1a74ccf5802 --- /dev/null +++ b/CalibPPS/AlignmentRelative/src/AlignmentAlgorithm.cc @@ -0,0 +1,16 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "CalibPPS/AlignmentRelative/interface/AlignmentAlgorithm.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentTask.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +//---------------------------------------------------------------------------------------------------- + +AlignmentAlgorithm::AlignmentAlgorithm(const edm::ParameterSet& ps, AlignmentTask* _t) + : verbosity(ps.getUntrackedParameter("verbosity", 0)), + task(_t), + singularLimit(ps.getParameter("singularLimit")) {} diff --git a/CalibPPS/AlignmentRelative/src/AlignmentGeometry.cc b/CalibPPS/AlignmentRelative/src/AlignmentGeometry.cc new file mode 100644 index 0000000000000..88971fcfaf50c --- /dev/null +++ b/CalibPPS/AlignmentRelative/src/AlignmentGeometry.cc @@ -0,0 +1,46 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h" +#include "CalibPPS/AlignmentRelative/interface/Utilities.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h" + +using namespace std; + +//---------------------------------------------------------------------------------------------------- + +const DetGeometry &AlignmentGeometry::get(unsigned int id) const { + auto it = sensorGeometry.find(id); + if (it == sensorGeometry.end()) + throw cms::Exception("PPS") << "No geometry available for sensor " << id << "."; + + return it->second; +} + +//---------------------------------------------------------------------------------------------------- + +void AlignmentGeometry::insert(unsigned int id, const DetGeometry &g) { sensorGeometry[id] = g; } + +//---------------------------------------------------------------------------------------------------- + +void AlignmentGeometry::print() const { + for (const auto &it : sensorGeometry) { + printId(it.first); + + printf(" z = %+10.4f mm │ shift: x = %+7.3f mm, y = %+7.3f mm │ ", it.second.z, it.second.sx, it.second.sy); + + for (const auto &dit : it.second.directionData) { + printf("dir%u: %+.3f, %+.3f, %+.3f │ ", dit.first, dit.second.dx, dit.second.dy, dit.second.dz); + } + + if (CTPPSDetId(it.first).subdetId() == CTPPSDetId::sdTrackingStrip) + printf("%s", (it.second.isU) ? "U-det" : "V-det"); + + printf("\n"); + } +} diff --git a/CalibPPS/AlignmentRelative/src/AlignmentTask.cc b/CalibPPS/AlignmentRelative/src/AlignmentTask.cc new file mode 100644 index 0000000000000..85dbe488d302e --- /dev/null +++ b/CalibPPS/AlignmentRelative/src/AlignmentTask.cc @@ -0,0 +1,547 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "CalibPPS/AlignmentRelative/interface/AlignmentTask.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentConstraint.h" + +#include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h" +#include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h" + +#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h" + +#include + +//---------------------------------------------------------------------------------------------------- + +using namespace std; +using namespace edm; + +//---------------------------------------------------------------------------------------------------- + +AlignmentTask::AlignmentTask(const ParameterSet &ps) + : resolveShR(ps.getParameter("resolveShR")), + resolveShZ(ps.getParameter("resolveShZ")), + resolveRotZ(ps.getParameter("resolveRotZ")), + + oneRotZPerPot(ps.getParameter("oneRotZPerPot")), + useEqualMeanUMeanVRotZConstraints(ps.getParameter("useEqualMeanUMeanVRotZConstraints")), + + fixedDetectorsConstraints(ps.getParameterSet("fixedDetectorsConstraints")), + standardConstraints(ps.getParameterSet("standardConstraints")) { + if (resolveShR) { + quantityClasses.push_back(qcShR1); + quantityClasses.push_back(qcShR2); + } + + if (resolveShZ) { + quantityClasses.push_back(qcShZ); + } + + if (resolveRotZ) { + quantityClasses.push_back(qcRotZ); + } +} + +//---------------------------------------------------------------------------------------------------- + +void AlignmentTask::buildGeometry(const vector &rpDecIds, + const vector &excludedSensors, + const CTPPSGeometry *input, + double z0, + AlignmentGeometry &geometry) { + geometry.z0 = z0; + + // traverse full known geometry + for (auto it = input->beginSensor(); it != input->endSensor(); ++it) { + // skip excluded sensors + if (find(excludedSensors.begin(), excludedSensors.end(), it->first) != excludedSensors.end()) + continue; + + // is RP selected? + const CTPPSDetId detId(it->first); + const unsigned int rpDecId = 100 * detId.arm() + 10 * detId.station() + detId.rp(); + if (find(rpDecIds.begin(), rpDecIds.end(), rpDecId) == rpDecIds.end()) + continue; + + // extract geometry data + CTPPSGeometry::Vector c = input->localToGlobal(detId, CTPPSGeometry::Vector(0., 0., 0.)); + CTPPSGeometry::Vector d1 = input->localToGlobal(detId, CTPPSGeometry::Vector(1., 0., 0.)) - c; + CTPPSGeometry::Vector d2 = input->localToGlobal(detId, CTPPSGeometry::Vector(0., 1., 0.)) - c; + + // for strips: is it U plane? + bool isU = false; + if (detId.subdetId() == CTPPSDetId::sdTrackingStrip) { + TotemRPDetId stripDetId(it->first); + unsigned int rpNum = stripDetId.rp(); + unsigned int plNum = stripDetId.plane(); + isU = (plNum % 2 != 0); + if (rpNum == 2 || rpNum == 3) + isU = !isU; + } + + DetGeometry dg(c.z() - z0, c.x(), c.y(), isU); + dg.setDirection(1, d1.x(), d1.y(), d1.z()); + dg.setDirection(2, d2.x(), d2.y(), d2.z()); + geometry.insert(it->first, dg); + } +} + +//---------------------------------------------------------------------------------------------------- + +void AlignmentTask::buildIndexMaps() { + // remove old mapping + mapMeasurementIndeces.clear(); + mapQuantityIndeces.clear(); + + // loop over all classes + for (const auto &qcl : quantityClasses) { + // create entry for this class + mapMeasurementIndeces[qcl]; + + // loop over all sensors + unsigned int idxMeas = 0; + unsigned int idxQuan = 0; + for (const auto &git : geometry.getSensorMap()) { + const unsigned int detId = git.first; + const unsigned int subdetId = CTPPSDetId(git.first).subdetId(); + + // update measurement map + if (qcl == qcShR1) { + if (subdetId == CTPPSDetId::sdTimingDiamond) + mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++; + if (subdetId == CTPPSDetId::sdTrackingPixel) + mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++; + } + + if (qcl == qcShR2) { + if (subdetId == CTPPSDetId::sdTrackingStrip) + mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++; + if (subdetId == CTPPSDetId::sdTrackingPixel) + mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++; + } + + if (qcl == qcShZ) { + if (subdetId == CTPPSDetId::sdTrackingStrip) + mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++; + if (subdetId == CTPPSDetId::sdTrackingPixel) + mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++; + if (subdetId == CTPPSDetId::sdTrackingPixel) + mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++; + if (subdetId == CTPPSDetId::sdTimingDiamond) + mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++; + } + + if (qcl == qcRotZ) { + if (subdetId == CTPPSDetId::sdTrackingStrip) + mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++; + if (subdetId == CTPPSDetId::sdTrackingPixel) + mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++; + if (subdetId == CTPPSDetId::sdTrackingPixel) + mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++; + if (subdetId == CTPPSDetId::sdTimingDiamond) + mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++; + } + + // update quantity map + if (qcl == qcShR1) { + if (subdetId == CTPPSDetId::sdTimingDiamond) + mapQuantityIndeces[qcl][detId] = idxQuan++; + if (subdetId == CTPPSDetId::sdTrackingPixel) + mapQuantityIndeces[qcl][detId] = idxQuan++; + } + + if (qcl == qcShR2) { + if (subdetId == CTPPSDetId::sdTrackingStrip) + mapQuantityIndeces[qcl][detId] = idxQuan++; + if (subdetId == CTPPSDetId::sdTrackingPixel) + mapQuantityIndeces[qcl][detId] = idxQuan++; + } + + if (qcl == qcShZ) { + if (subdetId == CTPPSDetId::sdTrackingStrip) + mapQuantityIndeces[qcl][detId] = idxQuan++; + if (subdetId == CTPPSDetId::sdTimingDiamond) + mapQuantityIndeces[qcl][detId] = idxQuan++; + if (subdetId == CTPPSDetId::sdTrackingPixel) + mapQuantityIndeces[qcl][detId] = idxQuan++; + } + + if (qcl == qcRotZ) { + if (subdetId == CTPPSDetId::sdTrackingStrip) + mapQuantityIndeces[qcl][detId] = idxQuan++; + if (subdetId == CTPPSDetId::sdTimingDiamond) + mapQuantityIndeces[qcl][detId] = idxQuan++; + if (subdetId == CTPPSDetId::sdTrackingPixel) + mapQuantityIndeces[qcl][detId] = idxQuan++; + } + } + } +} + +//---------------------------------------------------------------------------------------------------- + +signed int AlignmentTask::getMeasurementIndex(QuantityClass cl, unsigned int detId, unsigned int dirIdx) const { + auto clit = mapMeasurementIndeces.find(cl); + if (clit == mapMeasurementIndeces.end()) + return -1; + + auto it = clit->second.find({detId, dirIdx}); + if (it == clit->second.end()) + return -1; + + return it->second; +} + +//---------------------------------------------------------------------------------------------------- + +signed int AlignmentTask::getQuantityIndex(QuantityClass cl, unsigned int detId) const { + auto clit = mapQuantityIndeces.find(cl); + if (clit == mapQuantityIndeces.end()) + return -1; + + auto it = clit->second.find(detId); + if (it == clit->second.end()) + return -1; + + return it->second; +} + +//---------------------------------------------------------------------------------------------------- + +string AlignmentTask::quantityClassTag(QuantityClass qc) const { + switch (qc) { + case qcShR1: + return "ShR1"; + case qcShR2: + return "ShR2"; + case qcShZ: + return "ShZ"; + case qcRotZ: + return "RotZ"; + } + + throw cms::Exception("PPS") << "Unknown quantity class " << qc << "."; +} + +//---------------------------------------------------------------------------------------------------- + +unsigned int AlignmentTask::measurementsOfClass(QuantityClass qc) const { + auto it = mapMeasurementIndeces.find(qc); + if (it == mapMeasurementIndeces.end()) + return 0; + else + return it->second.size(); +} + +//---------------------------------------------------------------------------------------------------- + +unsigned int AlignmentTask::quantitiesOfClass(QuantityClass qc) const { + auto it = mapQuantityIndeces.find(qc); + if (it == mapQuantityIndeces.end()) + return 0; + else + return it->second.size(); +} + +//---------------------------------------------------------------------------------------------------- + +void AlignmentTask::buildFixedDetectorsConstraints(vector &constraints) const { + for (auto &quantityClass : quantityClasses) { + // get input + const string &tag = quantityClassTag(quantityClass); + + const ParameterSet &classSettings = fixedDetectorsConstraints.getParameterSet(tag.c_str()); + vector ids(classSettings.getParameter>("ids")); + vector values(classSettings.getParameter>("values")); + + if (ids.size() != values.size()) + throw cms::Exception("PPS") << "Different number of constraint ids and values for " << tag << "."; + + // determine number of constraints + unsigned int size = ids.size(); + + // just one basic constraint + if (oneRotZPerPot && quantityClass == qcRotZ) { + if (size > 1) + size = 1; + } + + // build constraints + for (unsigned int j = 0; j < size; j++) { + // prepare empty constraint + AlignmentConstraint ac; + + for (auto &qcit : quantityClasses) { + ac.coef[qcit].ResizeTo(quantitiesOfClass(qcit)); + ac.coef[qcit].Zero(); + } + + // set constraint name + char buf[40]; + sprintf(buf, "%s: fixed plane %4u", tag.c_str(), ids[j]); + ac.name = buf; + + // get quantity index + signed int qIndex = getQuantityIndex(quantityClass, ids[j]); + if (qIndex < 0) + throw cms::Exception("AlignmentTask::BuildFixedDetectorsConstraints") + << "Quantity index for class " << quantityClass << " and id " << ids[j] << " is " << qIndex; + + // set constraint coefficient and value + ac.coef[quantityClass][qIndex] = 1.; + ac.val = values[j] * 1E-3; + + // save constraint + constraints.push_back(ac); + } + + if (oneRotZPerPot && quantityClass == qcRotZ) + buildOneRotZPerPotConstraints(constraints); + } +} + +//---------------------------------------------------------------------------------------------------- + +void AlignmentTask::buildStandardConstraints(vector &constraints) const { + const vector &decUnitIds = standardConstraints.getParameter>("units"); + + // count planes in RPs + map planesPerPot; + for (const auto &it : geometry.getSensorMap()) { + CTPPSDetId detId(it.first); + planesPerPot[detId.rpId()]++; + } + + // ShR constraints + if (resolveShR) { + for (const auto &decUnitId : decUnitIds) { + // prepare empty constraints + AlignmentConstraint ac_X; + for (auto &qcit : quantityClasses) { + ac_X.coef[qcit].ResizeTo(quantitiesOfClass(qcit)); + ac_X.coef[qcit].Zero(); + } + ac_X.val = 0; + + AlignmentConstraint ac_Y(ac_X); + + // set constraint names + char buf[50]; + sprintf(buf, "ShR: unit %u, MeanX=0", decUnitId); + ac_X.name = buf; + sprintf(buf, "ShR: unit %u, MeanY=0", decUnitId); + ac_Y.name = buf; + + // traverse geometry + for (const auto &git : geometry.getSensorMap()) { + // stop is sensor not in the selected arm + CTPPSDetId senId(git.first); + unsigned int senDecUnit = senId.arm() * 100 + senId.station() * 10; + if (senId.rp() > 2) + senDecUnit += 1; + + if (senDecUnit != decUnitId) + continue; + + // fill constraint for strip sensors + if (senId.subdetId() == CTPPSDetId::sdTrackingStrip) { + signed int qIndex = getQuantityIndex(qcShR2, git.first); + if (qIndex < 0) + throw cms::Exception("AlignmentTask::BuildStandardConstraints") + << "Cannot get quantity index for class " << qcShR2 << " and sensor id " << git.first << "."; + + // determine weight + const double weight = 1. / planesPerPot[senId.rpId()]; + + // set constraint coefficients + ac_X.coef[qcShR2][qIndex] = git.second.getDirectionData(2).dx * weight; + ac_Y.coef[qcShR2][qIndex] = git.second.getDirectionData(2).dy * weight; + } + + // fill constraint for pixel sensors + if (senId.subdetId() == CTPPSDetId::sdTrackingPixel) { + // get quantity indeces + const signed int qIndex1 = getQuantityIndex(qcShR1, git.first); + if (qIndex1 < 0) + throw cms::Exception("AlignmentTask::BuildStandardConstraints") + << "Cannot get quantity index for class " << qcShR1 << " and sensor id " << git.first << "."; + + const signed int qIndex2 = getQuantityIndex(qcShR2, git.first); + if (qIndex2 < 0) + throw cms::Exception("AlignmentTask::BuildStandardConstraints") + << "Cannot get quantity index for class " << qcShR2 << " and sensor id " << git.first << "."; + + // determine weight (two constraints per plane) + const double weight = 0.5 / planesPerPot[senId.rpId()]; + + // get geometry + const double d1x = git.second.getDirectionData(1).dx; + const double d1y = git.second.getDirectionData(1).dy; + const double d2x = git.second.getDirectionData(2).dx; + const double d2y = git.second.getDirectionData(2).dy; + + // calculate coefficients, by inversion of this matrix relation + // [ s1 ] = [ d1x d1y ] * [ de x ] + // [ s2 ] [ d2x d2y ] [ de y ] + const double D = d1x * d2y - d1y * d2x; + const double coef_x_s1 = +d2y / D; + const double coef_y_s1 = -d2x / D; + const double coef_x_s2 = -d1y / D; + const double coef_y_s2 = +d1x / D; + + // set constraint coefficients + ac_X.coef[qcShR1][qIndex1] = coef_x_s1 * weight; + ac_Y.coef[qcShR1][qIndex1] = coef_y_s1 * weight; + ac_X.coef[qcShR2][qIndex2] = coef_x_s2 * weight; + ac_Y.coef[qcShR2][qIndex2] = coef_y_s2 * weight; + } + } + + // add constraints + constraints.push_back(ac_X); + constraints.push_back(ac_Y); + } + } + + // RotZ constraints + if (resolveRotZ) { + for (const auto &decUnitId : decUnitIds) { + // prepare empty constraints + AlignmentConstraint ac; + for (unsigned int i = 0; i < quantityClasses.size(); i++) { + ac.coef[quantityClasses[i]].ResizeTo(quantitiesOfClass(quantityClasses[i])); + ac.coef[quantityClasses[i]].Zero(); + } + ac.val = 0; + + char buf[50]; + sprintf(buf, "RotZ: unit %u, Mean=0", decUnitId); + ac.name = buf; + + // traverse geometry + for (const auto &git : geometry.getSensorMap()) { + // stop is sensor not in the selected arm + CTPPSDetId senId(git.first); + unsigned int senDecUnit = senId.arm() * 100 + senId.station() * 10; + if (senId.rp() > 2) + senDecUnit += 1; + + if (senDecUnit != decUnitId) + continue; + + // determine weight + const double weight = 1. / planesPerPot[senId.rpId()]; + + // set coefficient + signed int qIndex = getQuantityIndex(qcRotZ, git.first); + ac.coef[qcRotZ][qIndex] = weight; + } + + constraints.push_back(ac); + } + } + + if (resolveRotZ && oneRotZPerPot) + buildOneRotZPerPotConstraints(constraints); + + if (resolveRotZ && useEqualMeanUMeanVRotZConstraints) + buildEqualMeanUMeanVRotZConstraints(constraints); +} + +//---------------------------------------------------------------------------------------------------- + +void AlignmentTask::buildOneRotZPerPotConstraints(std::vector &constraints) const { + // build map rp id --> sensor ids + map> m; + for (const auto &p : geometry.getSensorMap()) { + CTPPSDetId detId(p.first); + CTPPSDetId rpId = detId.rpId(); + unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); + m[decRPId].push_back(p.first); + } + + // traverse all RPs + for (const auto &p : m) { + // build n_planes-1 constraints + unsigned int prev_detId = 0; + for (const auto &detId : p.second) { + if (prev_detId != 0) { + AlignmentConstraint ac; + + char buf[100]; + sprintf(buf, "RotZ: RP %u, plane %u = plane %u", p.first, prev_detId, detId); + ac.name = buf; + + ac.val = 0; + + for (auto &qcit : quantityClasses) { + ac.coef[qcit].ResizeTo(quantitiesOfClass(qcit)); + ac.coef[qcit].Zero(); + } + + signed int qIdx1 = getQuantityIndex(qcRotZ, prev_detId); + signed int qIdx2 = getQuantityIndex(qcRotZ, detId); + + ac.coef[qcRotZ][qIdx1] = +1.; + ac.coef[qcRotZ][qIdx2] = -1.; + + constraints.push_back(ac); + } + + prev_detId = detId; + } + } +} + +//---------------------------------------------------------------------------------------------------- + +void AlignmentTask::buildEqualMeanUMeanVRotZConstraints(vector &constraints) const { + // build map strip rp id --> pair( vector of U planes, vector of V planes ) + map, vector>> m; + for (const auto &p : geometry.getSensorMap()) { + CTPPSDetId detId(p.first); + + if (detId.subdetId() != CTPPSDetId::sdTrackingStrip) + continue; + + CTPPSDetId rpId = detId.rpId(); + unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); + + if (p.second.isU) + m[decRPId].first.push_back(p.first); + else + m[decRPId].second.push_back(p.first); + } + + // loop over RPs + for (const auto &p : m) { + AlignmentConstraint ac; + + char buf[100]; + sprintf(buf, "RotZ: RP %u, MeanU = MeanV", p.first); + ac.name = buf; + + ac.val = 0; + + for (auto &qcit : quantityClasses) { + ac.coef[qcit].ResizeTo(quantitiesOfClass(qcit)); + ac.coef[qcit].Zero(); + } + + for (const string &proj : {"U", "V"}) { + const auto &planes = (proj == "U") ? p.second.first : p.second.second; + const double c = ((proj == "U") ? -1. : +1.) / planes.size(); + + for (const auto &plane : planes) { + signed int qIdx = getQuantityIndex(qcRotZ, plane); + ac.coef[qcRotZ][qIdx] = c; + + TotemRPDetId plId(plane); + } + } + + constraints.push_back(ac); + } +} diff --git a/CalibPPS/AlignmentRelative/src/IdealResult.cc b/CalibPPS/AlignmentRelative/src/IdealResult.cc new file mode 100644 index 0000000000000..b7cf365f7e9ee --- /dev/null +++ b/CalibPPS/AlignmentRelative/src/IdealResult.cc @@ -0,0 +1,302 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "CalibPPS/AlignmentRelative/interface/IdealResult.h" + +#include "CalibPPS/AlignmentRelative/interface/AlignmentTask.h" + +#include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h" +#include "Geometry/Records/interface/VeryForwardMisalignedGeometryRecord.h" + +#include "TMatrixD.h" +#include "TVectorD.h" + +using namespace std; + +//---------------------------------------------------------------------------------------------------- + +IdealResult::IdealResult(const edm::ParameterSet &gps, AlignmentTask *_t) : AlignmentAlgorithm(gps, _t) {} + +//---------------------------------------------------------------------------------------------------- + +void IdealResult::begin(const CTPPSGeometry *geometryReal, const CTPPSGeometry *geometryMisaligned) { + gReal = geometryReal; + gMisaligned = geometryMisaligned; +} +//---------------------------------------------------------------------------------------------------- + +unsigned int IdealResult::solve(const std::vector &constraints, + std::map &results, + TDirectory *dir) { + /* + 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 offsets; + map 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 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; +} diff --git a/CalibPPS/AlignmentRelative/src/JanAlignmentAlgorithm.cc b/CalibPPS/AlignmentRelative/src/JanAlignmentAlgorithm.cc new file mode 100644 index 0000000000000..4d38f8b2e8a49 --- /dev/null +++ b/CalibPPS/AlignmentRelative/src/JanAlignmentAlgorithm.cc @@ -0,0 +1,772 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "CalibPPS/AlignmentRelative/interface/JanAlignmentAlgorithm.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentTask.h" +#include "CalibPPS/AlignmentRelative/interface/Utilities.h" + +#include "TMatrixDSymEigen.h" +#include "TDecompSVD.h" +#include "TFile.h" +#include "TCanvas.h" +#include "TH2D.h" + +#include + +//#define DEBUG 1 + +using namespace std; +using namespace edm; + +//---------------------------------------------------------------------------------------------------- + +JanAlignmentAlgorithm::JanAlignmentAlgorithm(const ParameterSet &ps, AlignmentTask *_t) + : AlignmentAlgorithm(ps, _t), Sc(nullptr), Mc(nullptr) { + const ParameterSet &lps = ps.getParameterSet("JanAlignmentAlgorithm"); + weakLimit = lps.getParameter("weakLimit"); + stopOnSingularModes = lps.getParameter("stopOnSingularModes"); + buildDiagnosticPlots = lps.getParameter("buildDiagnosticPlots"); +} + +//---------------------------------------------------------------------------------------------------- + +JanAlignmentAlgorithm::~JanAlignmentAlgorithm() {} + +//---------------------------------------------------------------------------------------------------- + +void JanAlignmentAlgorithm::begin(const CTPPSGeometry *geometryReal, const CTPPSGeometry *geometryMisaligned) { + // initialize M and S components + Mc = new TVectorD[task->quantityClasses.size()]; + Sc = new TMatrixD *[task->quantityClasses.size()]; + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + unsigned int rows = task->quantitiesOfClass(task->quantityClasses[i]); + + Mc[i].ResizeTo(rows); + Mc[i].Zero(); + + Sc[i] = new TMatrixD[task->quantityClasses.size()]; + for (unsigned int j = 0; j < task->quantityClasses.size(); j++) { + unsigned int cols = task->quantitiesOfClass(task->quantityClasses[j]); + Sc[i][j].ResizeTo(rows, cols); + Sc[i][j].Zero(); + } + } + + // prepare statistics plots + if (buildDiagnosticPlots) { + for (const auto &it : task->geometry.getSensorMap()) { + unsigned int id = it.first; + char buf[50]; + DetStat s; + + sprintf(buf, "%u: m distribution", id); + s.m_dist = new TH1D(buf, ";u or v (mm)", 100, -25., 25.); + + sprintf(buf, "%u: R distribution", id); + s.R_dist = new TH1D(buf, ";R (mm)", 500, -0.5, 0.5); + + for (unsigned int c = 0; c < task->quantityClasses.size(); c++) { + sprintf(buf, "%u: coef, %s", id, task->quantityClassTag(task->quantityClasses[c]).c_str()); + s.coefHist.push_back(new TH1D(buf, ";coefficient", 100, -2., +2.)); + + sprintf(buf, "%u: R vs. coef, %s", id, task->quantityClassTag(task->quantityClasses[c]).c_str()); + TGraph *g = new TGraph(); + g->SetName(buf); + g->SetTitle(";coefficient;residual (mm)"); + s.resVsCoef.push_back(g); + } + + statistics[id] = s; + } + } + + events = 0; +} + +//---------------------------------------------------------------------------------------------------- + +void JanAlignmentAlgorithm::feed(const HitCollection &selection, const LocalTrackFit &trackFit) { + if (verbosity > 9) + printf("\n>> JanAlignmentAlgorithm::Feed\n"); + + events++; + + // prepare fit - make z0 compatible + double hax = trackFit.ax; + double hay = trackFit.ay; + double hbx = trackFit.bx + trackFit.ax * (task->geometry.z0 - trackFit.z0); + double hby = trackFit.by + trackFit.ay * (task->geometry.z0 - trackFit.z0); + + // prepare Gamma matrices (full of zeros) + TMatrixD *Ga = new TMatrixD[task->quantityClasses.size()]; + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + Ga[i].ResizeTo(selection.size(), Mc[i].GetNrows()); + Ga[i].Zero(); + } + + TMatrixD A(selection.size(), 4); + TMatrixD Vi(selection.size(), selection.size()); + TVectorD m(selection.size()); + + set rpSet; + if (buildDiagnosticPlots) { + for (const auto &hit : selection) { + CTPPSDetId detId(hit.id); + const unsigned int rpDecId = 100 * detId.arm() + 10 * detId.station() + detId.rp(); + rpSet.insert(rpDecId); + } + } + + // fill fit matrix and Gamma matrices + unsigned int j = 0; + + for (HitCollection::const_iterator hit = selection.begin(); hit != selection.end(); ++hit, ++j) { + unsigned int id = hit->id; + + const DetGeometry &d = task->geometry.get(id); + const auto &dirData = d.getDirectionData(hit->dirIdx); + + A(j, 0) = hit->z * dirData.dx; + A(j, 1) = dirData.dx; + A(j, 2) = hit->z * dirData.dy; + A(j, 3) = dirData.dy; + + m(j) = hit->position + dirData.s - (hit->z - d.z) * dirData.dz; // in mm + + Vi(j, j) = 1. / hit->sigma / hit->sigma; + + double C = dirData.dx, S = dirData.dy; + + double hx = hax * hit->z + hbx; // in mm + double hy = hay * hit->z + hby; + double R = m(j) - (hx * C + hy * S); // (standard) residual + + if (buildDiagnosticPlots) { + statistics[id].m_dist->Fill(m(j)); + statistics[id].R_dist->Fill(R); + } + + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + // check compatibility + signed int matrixIndex = task->getMeasurementIndex(task->quantityClasses[i], hit->id, hit->dirIdx); + if (matrixIndex < 0) + continue; + + matrixIndex = task->getQuantityIndex(task->quantityClasses[i], hit->id); + + switch (task->quantityClasses[i]) { + case AlignmentTask::qcShR1: + Ga[i][j][matrixIndex] = -1.; + break; + + case AlignmentTask::qcShR2: + Ga[i][j][matrixIndex] = -1.; + break; + + case AlignmentTask::qcShZ: + Ga[i][j][matrixIndex] = hax * C + hay * S; + break; + + case AlignmentTask::qcRotZ: + Ga[i][j][matrixIndex] = (hax * hit->z + hbx - d.sx) * (-S) + (hay * hit->z + hby - d.sy) * C; + break; + } + + if (buildDiagnosticPlots) { + double c = Ga[i][j][matrixIndex]; + DetStat &s = statistics[id]; + s.coefHist[i]->Fill(c); + s.resVsCoef[i]->SetPoint(s.resVsCoef[i]->GetN(), c, R); + + if (task->quantityClasses[i] == AlignmentTask::qcRotZ) { + map, ScatterPlot>::iterator hit = s.resVsCoefRot_perRPSet.find(rpSet); + if (hit == s.resVsCoefRot_perRPSet.end()) { + ScatterPlot sp; + sp.g = new TGraph(); + sp.h = new TH2D("", "", 40, -20., +20., 60, -0.15, +0.15); + hit = s.resVsCoefRot_perRPSet.insert(pair, ScatterPlot>(rpSet, sp)).first; + } + hit->second.g->SetPoint(hit->second.g->GetN(), c, R); + hit->second.h->Fill(c, R); + } + } + } + } + + // sigma matrix + TMatrixD AT(TMatrixD::kTransposed, A); + TMatrixD ATViA(4, 4); + ATViA = AT * Vi * A; + TMatrixD ATViAI(ATViA); + ATViAI = ATViA.Invert(); + TMatrixD sigma(Vi); + sigma -= Vi * A * ATViAI * AT * Vi; + + // traspose Gamma matrices + TMatrixD *GaT = new TMatrixD[task->quantityClasses.size()]; + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + GaT[i].ResizeTo(Mc[i].GetNrows(), selection.size()); + GaT[i].Transpose(Ga[i]); + } + + // normalized residuals + TVectorD r(selection.size()); + r = sigma * m; + + // increment M + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + if (Mc[i].GetNrows() < 1) + continue; + + Mc[i] += GaT[i] * r; + } + + // increment S + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + for (unsigned int j = 0; j < task->quantityClasses.size(); j++) { + if (Sc[i][j].GetNrows() < 1 || Sc[i][j].GetNcols() < 1) + continue; + + Sc[i][j] += GaT[i] * sigma * Ga[j]; + } + } + +#ifdef DEBUG + printf("* checking normalized residuals, selection.size = %u\n", selection.size()); + r.Print(); + + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + printf("- class %u\n", i); + TVectorD t(Mc[i].GetNrows()); + for (int j = 0; j < t.GetNrows(); j++) + t[j] = 1.; + t.Print(); + + Ga[i].Print(); + + TVectorD tt(selection.size()); + tt = sigma * Ga[i] * t; + + double ttn = sqrt(tt.Norm2Sqr()); + printf("|tt| = %E\n", ttn); + if (ttn > 1E-8) + tt.Print(); + } +#endif + + delete[] Ga; + delete[] GaT; +} + +//---------------------------------------------------------------------------------------------------- + +void JanAlignmentAlgorithm::analyze() { + if (verbosity > 2) + printf("\n>> JanAlignmentAlgorithm::Analyze\n"); + + // calculate full dimension + unsigned int dim = 0; + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) + dim += Mc[i].GetNrows(); + + if (verbosity > 2) { + printf("\tsensors: %u\n", task->geometry.getNumberOfDetectors()); + printf("\tfull dimension: %u\n", dim); + printf("\tquantity classes: %lu\n", task->quantityClasses.size()); + } + + // build full M + M.ResizeTo(dim); + unsigned int offset = 0; + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + M.SetSub(offset, Mc[i]); + offset += Mc[i].GetNrows(); + } + + // build full S + S.ResizeTo(dim, dim); + unsigned int r_offset = 0, c_offset = 0; + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + c_offset = 0; + unsigned int r_size = 0, c_size = 0; + for (unsigned int j = 0; j < task->quantityClasses.size(); j++) { + r_size = Sc[i][j].GetNrows(); + c_size = Sc[i][j].GetNcols(); + + if (r_size < 1 || c_size < 1) + continue; + + TMatrixDSub(S, r_offset, r_offset + r_size - 1, c_offset, c_offset + c_size - 1) = Sc[i][j]; + c_offset += c_size; + } + r_offset += r_size; + } + + // analyze symmetricity + if (verbosity >= 3) { + double maxDiff = 0., maxElem = 0.; + for (unsigned int i = 0; i < dim; i++) { + for (unsigned int j = 0; j < dim; j++) { + double diff = S[i][j] - S[j][i]; + if (fabs(diff) > maxDiff) + maxDiff = diff; + if (S[i][j] > maxElem) + maxElem = S[i][j]; + } + } + + printf("\n* S matrix:\n\tdimension = %i\n\tmaximum asymmetry: %E\t(ratio to maximum element %E)\n", + dim, + maxDiff, + maxDiff / maxElem); + } + + // make a symmetric copy + TMatrixDSym S_sym(dim); + for (unsigned int j = 0; j < dim; j++) { + for (unsigned int i = 0; i < dim; i++) { + S_sym[i][j] = S[i][j]; + } + } + + // eigen analysis of S + TMatrixDSymEigen S_eig(S_sym); + const TVectorD &S_eigVal_temp = S_eig.GetEigenValues(); + S_eigVal.ResizeTo(S_eigVal_temp.GetNrows()); + S_eigVal = S_eigVal_temp; + const TMatrixD &S_eigVec_temp = S_eig.GetEigenVectors(); + S_eigVec.ResizeTo(S_eigVec_temp); + S_eigVec = S_eigVec_temp; + + // identify singular modes + for (int i = 0; i < S_eigVal.GetNrows(); i++) { + double nev = S_eigVal[i] / events; + if (fabs(nev) < singularLimit) { + SingularMode sM{S_eigVal[i], TMatrixDColumn(S_eigVec, i), i}; + singularModes.push_back(sM); + } + } + +#if DEBUG + // print singular vectors + if (singularModes.size() > 0) { + printf("\n* S singular modes\n | "); + for (unsigned int i = 0; i < singularModes.size(); i++) + printf("%+10.3E ", singularModes[i].val); + printf("\n-- | "); + + for (unsigned int i = 0; i < singularModes.size(); i++) + printf("---------- "); + printf("\n"); + + for (unsigned int j = 0; j < dim; j++) { + printf("%2u | ", j); + for (unsigned int i = 0; i < singularModes.size(); i++) { + printf("%+10.3E ", singularModes[i].vec[j]); + } + printf("\n"); + } + } else + printf("\n* S has no singular modes\n"); +#endif +} + +//---------------------------------------------------------------------------------------------------- + +unsigned int JanAlignmentAlgorithm::solve(const std::vector &constraints, + map &results, + TDirectory *dir) { + if (verbosity) + printf(">> JanAlignmentAlgorithm::Solve\n"); + + results.clear(); + + // build C matrix + unsigned int dim = S.GetNrows(); + TMatrixD C(dim, constraints.size()); + TMatrixD C2(dim, constraints.size()); + for (unsigned int i = 0; i < constraints.size(); i++) { + 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] = events * cv[k]; + C2[offset][i] = events * cv[k] * 1E3; + offset++; + } + } + } + +#ifdef DEBUG + printf("\n* constraint matrix\n"); + Print(C); +#endif + + // build E matrix (singular vectors of S as its columns) + TMatrixD E(S.GetNrows(), singularModes.size()); + for (unsigned int i = 0; i < singularModes.size(); i++) + for (int j = 0; j < S.GetNrows(); j++) + E(j, i) = singularModes[i].vec[j]; + + // build CS matrix + TMatrixDSym CS(dim + constraints.size()); + TMatrixDSym CS2(dim + constraints.size()); + CS.Zero(); + CS2.Zero(); + for (unsigned int j = 0; j < dim; j++) { + for (unsigned int i = 0; i < dim; i++) { + CS[i][j] = S[i][j]; + CS2[i][j] = S[i][j]; + } + } + + for (unsigned int i = 0; i < constraints.size(); i++) { + for (unsigned int j = 0; j < dim; j++) { + CS[j][dim + i] = CS[dim + i][j] = C(j, i); + CS2[j][dim + i] = CS2[dim + i][j] = C2(j, i); + } + } + + // eigen analysis of CS matrix + TMatrixDSymEigen CS_eig(CS); + TVectorD CS_eigVal = CS_eig.GetEigenValues(); + TMatrixD CS_eigVec = CS_eig.GetEigenVectors(); + + // check regularity of CS matrix + if (verbosity >= 2) { + printf("\n* eigen values of CS and S matrices (events = %u)\n", events); + printf(" # CS norm. CS S norm. S\n"); + } + + unsigned int singularModeCount = 0; + vector weakModeIdx; + for (int i = 0; i < CS_eigVal.GetNrows(); i++) { + const double CS_nev = CS_eigVal[i] / events; + + if (fabs(CS_nev) < singularLimit) + singularModeCount++; + + if (verbosity >= 2) { + printf("%4i%+12.2E%+12.2E", i, CS_eigVal[i], CS_nev); + if (fabs(CS_nev) < singularLimit) { + singularModeCount++; + printf(" (S)"); + } else { + if (fabs(CS_nev) < weakLimit) { + weakModeIdx.push_back(i); + printf(" (W)"); + } else { + printf(" "); + } + } + + if (i < S_eigVal.GetNrows()) { + double S_nev = S_eigVal[i] / events; + printf("%+12.2E%+12.2E", S_eigVal[i], S_nev); + if (fabs(S_nev) < singularLimit) + printf(" (S)"); + else if (fabs(S_nev) < weakLimit) + printf(" (W)"); + } + + printf("\n"); + } + } + + if (verbosity >= 2) { + // print weak vectors + if (!weakModeIdx.empty()) { + unsigned int columns = 10; + unsigned int first = 0; + + while (first < weakModeIdx.size()) { + unsigned int last = first + columns; + if (last >= weakModeIdx.size()) + last = weakModeIdx.size(); + + printf("\n* CS weak modes\n | "); + for (unsigned int i = first; i < last; i++) + printf("%+10.3E ", CS_eigVal[weakModeIdx[i]]); + printf("\n--- | "); + + for (unsigned int i = first; i < last; i++) + printf("---------- "); + printf("\n"); + + // determine maximum elements + vector maxs; + for (unsigned int i = first; i < last; i++) { + double max = 0; + for (unsigned int j = 0; j < dim + constraints.size(); j++) { + double v = fabs(CS_eigVec(weakModeIdx[i], j)); + if (v > max) + max = v; + } + maxs.push_back(max); + } + + for (unsigned int j = 0; j < dim + constraints.size(); j++) { + printf("%3u | ", j); + for (unsigned int i = first; i < last; i++) { + double v = CS_eigVec(weakModeIdx[i], j); + if (fabs(v) / maxs[i - first] > 1E-3) + printf("%+10.3E ", v); + else + printf(" . "); + } + printf("\n"); + } + + first = last; + } + } else + printf("\n* CS has no weak modes\n"); + } + + // check the regularity of C^T E + if (verbosity >= 2) { + if (E.GetNcols() == C.GetNcols()) { + TMatrixD CTE(C, TMatrixD::kTransposeMult, E); + print(CTE, "* CTE matrix:"); + const double &det = CTE.Determinant(); + printf( + "\n* det(CTE) = %E, max(CTE) = %E, det(CTE)/max(CTE) = %E\n\tmax(C) = %E, max(E) = %E, " + "det(CTE)/max(C)/max(E) = %E\n", + det, + CTE.Max(), + det / CTE.Max(), + C.Max(), + E.Max(), + det / C.Max() / E.Max()); + } else { + printf(">> JanAlignmentAlgorithm::Solve > WARNING: C matrix has %u, while E matrix %u columns.\n", + C.GetNcols(), + E.GetNcols()); + } + } + + // stop if CS is singular + if (singularModeCount > 0 && stopOnSingularModes) { + LogError("PPS") << "\n>> JanAlignmentAlgorithm::Solve > ERROR: There are " << singularModeCount + << " singular modes in CS matrix. Stopping."; + return 1; + } + + // build MV vector + TVectorD MV(dim + constraints.size()); + for (unsigned int i = 0; i < dim; i++) + MV[i] = M[i]; + for (unsigned int i = 0; i < constraints.size(); i++) + MV[dim + i] = events * constraints[i].val; + + // perform inversion and solution + TMatrixD CSI(TMatrixD::kInverted, CS); + TMatrixD CS2I(TMatrixD::kInverted, CS2); + TVectorD AL(MV); + AL = CSI * MV; + + // evaluate error matrix + TMatrixD S0(S); // new parts full of zeros + S0.ResizeTo(dim + constraints.size(), dim + constraints.size()); + TMatrixD EM(CS); + EM = CSI * S0 * CSI; + + TMatrixD EM2(CS2); + EM2 = CS2I * S0 * CS2I; + + TMatrixD EMdiff(EM2 - EM); + + if (verbosity >= 3) { + double max1 = -1., max2 = -1., maxDiff = -1.; + for (int i = 0; i < EMdiff.GetNrows(); i++) { + for (int j = 0; j < EMdiff.GetNcols(); j++) { + if (maxDiff < fabs(EMdiff(i, j))) + maxDiff = fabs(EMdiff(i, j)); + + if (max1 < fabs(EM(i, j))) + max1 = fabs(EM(i, j)); + + if (max2 < fabs(EM2(i, j))) + max2 = fabs(EM2(i, j)); + } + } + + printf("EM max = %E, EM2 max = %E, EM2 - EM max = %E\n", max1, max2, maxDiff); + } + + // tests + TMatrixD &U = CS_eigVec; + TMatrixD UT(TMatrixD::kTransposed, U); + //TMatrixD CSEi(CS); + //CSEi = UT * CS * U; + //Print(CSEi, "CSEi"); + + TMatrixD EMEi(EM); + EMEi = UT * EM * U; + //Print(EMEi, "*EMEi"); + + if (verbosity >= 3) { + double max = -1.; + for (int i = 0; i < EMEi.GetNrows(); i++) { + for (int j = 0; j < EMEi.GetNcols(); j++) { + if (max < EMEi(i, j)) + max = EMEi(i, j); + } + } + + printf("max = %E\n", max); + } + + // print lambda values + if (verbosity >= 3) { + printf("\n* Lambda (from the contribution of singular modes to MV)\n"); + for (unsigned int i = 0; i < constraints.size(); i++) { + printf("\t%u (%25s)\t%+10.1E +- %10.1E\n", + i, + constraints[i].name.c_str(), + AL[dim + i] * 1E3, + sqrt(EM[i + dim][i + dim]) * 1E3); + } + } + + // fill results + unsigned int offset = 0; + vector offsets; + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + offsets.push_back(offset); + offset += Mc[i].GetNrows(); + } + + for (const auto &dit : task->geometry.getSensorMap()) { + AlignmentResult r; + + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + signed idx = task->getQuantityIndex(task->quantityClasses[i], dit.first); + if (idx < 0) + continue; + + unsigned int fi = offsets[i] + idx; + double v = AL[fi]; + double e = sqrt(EM[fi][fi]); + switch (task->quantityClasses[i]) { + case AlignmentTask::qcShR1: + r.setShR1(v); + r.setShR1Unc(e); + break; + case AlignmentTask::qcShR2: + r.setShR2(v); + r.setShR2Unc(e); + break; + case AlignmentTask::qcShZ: + r.setShZ(v); + r.setShZUnc(e); + break; + case AlignmentTask::qcRotZ: + r.setRotZ(v); + r.setRotZUnc(e); + break; + } + } + + results[dit.first] = r; + } + + // save matrices, eigen data, ... + if (dir) { + dir->cd(); + + S.Write("S"); + S_eigVal.Write("S_eigen_values"); + S_eigVec.Write("S_eigen_vectors"); + + E.Write("E"); + C.Write("C"); + + CS.Write("CS"); + CS_eigVal.Write("CS_eigen_values"); + CS_eigVec.Write("CS_eigen_vectors"); + + MV.Write("MV"); + AL.Write("AL"); + + S0.Write("S0"); + EM.Write("EM"); + } + + // success + return 0; +} + +//---------------------------------------------------------------------------------------------------- + +void JanAlignmentAlgorithm::end() { + delete[] Mc; + + for (unsigned int i = 0; i < task->quantityClasses.size(); i++) { + delete[] Sc[i]; + } + delete[] Sc; +} + +//---------------------------------------------------------------------------------------------------- + +void JanAlignmentAlgorithm::saveDiagnostics(TDirectory *dir) { + if (!buildDiagnosticPlots) + return; + + for (map::iterator it = statistics.begin(); it != statistics.end(); ++it) { + char buf[50]; + sprintf(buf, "%u", it->first); + gDirectory = dir->mkdir(buf); + + it->second.m_dist->Write(); + it->second.R_dist->Write(); + + for (unsigned int c = 0; c < task->quantityClasses.size(); c++) { + it->second.coefHist[c]->Write(); + it->second.resVsCoef[c]->Write(); + } + + gDirectory = gDirectory->mkdir("R vs. rot. coef, per RP set"); + TCanvas *c = new TCanvas; + c->SetName("R vs. rot. coef, overlapped"); + TH2D *frame = new TH2D("frame", "frame", 100, -20., +20., 100, -0.15, +0.15); + frame->Draw(); + unsigned int idx = 0; + for (map, ScatterPlot>::iterator iit = it->second.resVsCoefRot_perRPSet.begin(); + iit != it->second.resVsCoefRot_perRPSet.end(); + ++iit, ++idx) { + string label; + bool first = true; + for (set::iterator sit = iit->first.begin(); sit != iit->first.end(); ++sit) { + char buf[50]; + sprintf(buf, "%u", *sit); + + if (first) { + label = buf; + first = false; + } else { + label = label + ", " + buf; + } + } + + iit->second.g->SetTitle(";rotation coefficient (mm);residual (mm)"); + iit->second.g->SetMarkerColor(idx + 1); + iit->second.g->SetName(label.c_str()); + iit->second.g->Draw((idx == 0) ? "p" : "p"); + iit->second.g->Write(); + + iit->second.h->SetName((label + " (hist)").c_str()); + iit->second.h->SetTitle(";rotation coefficient (mm);residual (mm)"); + iit->second.h->Write(); + } + + gDirectory->cd(".."); + c->Write(); + } +} diff --git a/CalibPPS/AlignmentRelative/src/LocalTrackFitter.cc b/CalibPPS/AlignmentRelative/src/LocalTrackFitter.cc new file mode 100644 index 0000000000000..dc7042dfceb53 --- /dev/null +++ b/CalibPPS/AlignmentRelative/src/LocalTrackFitter.cc @@ -0,0 +1,223 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "CalibPPS/AlignmentRelative/interface/LocalTrackFitter.h" + +#include "CalibPPS/AlignmentRelative/interface/Utilities.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h" +#include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h" + +#include + +#include "TMatrixD.h" +#include "TVectorD.h" + +using namespace std; + +//---------------------------------------------------------------------------------------------------- + +LocalTrackFitter::LocalTrackFitter(const edm::ParameterSet &ps) + : verbosity(ps.getUntrackedParameter("verbosity", 0)), + minimumHitsPerProjectionPerRP(ps.getParameter("minimumHitsPerProjectionPerRP")), + maxResidualToSigma(ps.getParameter("maxResidualToSigma")) {} + +//---------------------------------------------------------------------------------------------------- + +bool LocalTrackFitter::fit(HitCollection &selection, const AlignmentGeometry &geometry, LocalTrackFit &trackFit) const { + if (verbosity > 5) + printf(">> LocalTrackFitter::Fit\n"); + + bool selectionChanged = true; + unsigned int loopCounter = 0; + while (selectionChanged) { + // fit/outlier-removal loop + while (selectionChanged) { + if (verbosity > 5) + printf("* fit loop %u\n", loopCounter++); + + bool fitFailed = false; + fitAndRemoveOutliers(selection, geometry, trackFit, fitFailed, selectionChanged); + + if (fitFailed) { + if (verbosity > 5) + printf("\tFIT FAILED\n"); + return false; + } + } + + // remove pots with too few active planes + removeInsufficientPots(selection, selectionChanged); + } + + return true; +} + +//---------------------------------------------------------------------------------------------------- + +void LocalTrackFitter::fitAndRemoveOutliers(HitCollection &selection, + const AlignmentGeometry &geometry, + LocalTrackFit &trackFit, + bool &failed, + bool &selectionChanged) const { + if (verbosity > 5) + printf(" - LocalTrackFitter::FitAndRemoveOutliers\n"); + + if (selection.empty()) { + failed = true; + return; + } + + // check if input size is sufficient + if (selection.size() < 4) { + failed = true; + return; + } + + // build matrices and vectors + TMatrixD A(selection.size(), 4); + TMatrixD Vi(selection.size(), selection.size()); + TVectorD measVec(selection.size()); + unsigned int j = 0; + for (auto it = selection.begin(); it != selection.end(); ++it, ++j) { + const unsigned int &detId = it->id; + + const DetGeometry &d = geometry.get(detId); + const auto &dirData = d.getDirectionData(it->dirIdx); + + A(j, 0) = it->z * dirData.dx; + A(j, 1) = dirData.dx; + A(j, 2) = it->z * dirData.dy; + A(j, 3) = dirData.dy; + + measVec(j) = it->position + dirData.s - (it->z - d.z) * dirData.dz; // in mm + + Vi(j, j) = 1. / it->sigma / it->sigma; + } + + // evaluate local track parameter estimates (h stands for hat) + TMatrixD AT(4, selection.size()); + AT.Transpose(A); + TMatrixD ATViA(4, 4); + ATViA = AT * Vi * A; + TMatrixD ATViAI(ATViA); + + try { + ATViAI = ATViA.Invert(); + } catch (cms::Exception &e) { + failed = true; + return; + } + + TVectorD theta(4); + theta = ATViAI * AT * Vi * measVec; + + // residuals + TVectorD R(measVec); + R -= A * theta; + + // save results to trackFit + trackFit.ax = theta(0); + trackFit.bx = theta(1); + trackFit.ay = theta(2); + trackFit.by = theta(3); + trackFit.z0 = geometry.z0; + trackFit.ndf = selection.size() - 4; + trackFit.chi_sq = 0; + for (int i = 0; i < R.GetNrows(); i++) + trackFit.chi_sq += R(i) * R(i) * Vi(i, i); + + if (verbosity > 5) { + printf(" ax = %.3f mrad, bx = %.4f mm, ay = %.3f mrad, by = %.4f mm, z0 = %.3f mm\n", + trackFit.ax * 1E3, + trackFit.bx, + trackFit.ay * 1E3, + trackFit.by, + trackFit.z0); + printf(" ndof = %i, chi^2/ndof/si^2 = %.3f\n", trackFit.ndf, trackFit.chi_sq / trackFit.ndf); + } + + // check residuals + selectionChanged = false; + TVectorD interpolation(A * theta); + j = 0; + for (auto it = selection.begin(); it != selection.end(); ++j) { + if (verbosity > 5) { + printf(" %2u, ", j); + printId(it->id); + printf(", dirIdx=%u: interpol = %+8.1f um, residual = %+6.1f um, residual / sigma = %+6.2f\n", + it->dirIdx, + interpolation[j] * 1E3, + R[j] * 1E3, + R[j] / it->sigma); + } + + double resToSigma = R[j] / it->sigma; + if (fabs(resToSigma) > maxResidualToSigma) { + selection.erase(it); + selectionChanged = true; + if (verbosity > 5) + printf(" Removed\n"); + } else + ++it; + } +} + +//---------------------------------------------------------------------------------------------------- + +void LocalTrackFitter::removeInsufficientPots(HitCollection &selection, bool &selectionChanged) const { + if (verbosity > 5) + printf(" - RemoveInsufficientPots\n"); + + selectionChanged = false; + + // map: RP id -> (active planes in projection 1, active planes in projection 2) + map, set > > planeMap; + for (auto it = selection.begin(); it != selection.end(); ++it) { + CTPPSDetId senId(it->id); + const unsigned int rpId = senId.rpId(); + + if (senId.subdetId() == CTPPSDetId::sdTrackingStrip) { + const unsigned int plane = TotemRPDetId(it->id).plane(); + if ((plane % 2) == 0) + planeMap[rpId].first.insert(senId); + else + planeMap[rpId].second.insert(senId); + } + + if (senId.subdetId() == CTPPSDetId::sdTrackingPixel) { + planeMap[rpId].first.insert(senId); + planeMap[rpId].second.insert(senId); + } + } + + // remove RPs with insufficient information + selectionChanged = false; + + for (auto it = planeMap.begin(); it != planeMap.end(); ++it) { + if (it->second.first.size() < minimumHitsPerProjectionPerRP || + it->second.second.size() < minimumHitsPerProjectionPerRP) { + if (verbosity > 5) + printf("\tRP %u: projection1 = %lu, projection 2 = %lu\n", + it->first, + it->second.first.size(), + it->second.second.size()); + + // remove all hits from that RP + for (auto dit = selection.begin(); dit != selection.end();) { + if (it->first == CTPPSDetId(dit->id).rpId()) { + if (verbosity > 5) + printf("\t\tremoving %u\n", dit->id); + selection.erase(dit); + selectionChanged = true; + } else { + dit++; + } + } + } + } +} diff --git a/CalibPPS/AlignmentRelative/src/StraightTrackAlignment.cc b/CalibPPS/AlignmentRelative/src/StraightTrackAlignment.cc new file mode 100644 index 0000000000000..55d94d8127778 --- /dev/null +++ b/CalibPPS/AlignmentRelative/src/StraightTrackAlignment.cc @@ -0,0 +1,1106 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsMethods.h" +#include "CondFormats/AlignmentRecord/interface/RPRealAlignmentRecord.h" +#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h" +#include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h" + +#include "CalibPPS/AlignmentRelative/interface/StraightTrackAlignment.h" +#include "CalibPPS/AlignmentRelative/interface/Utilities.h" + +#include "CalibPPS/AlignmentRelative/interface/IdealResult.h" +#include "CalibPPS/AlignmentRelative/interface/JanAlignmentAlgorithm.h" + +#include +#include +#include +#include + +#include "TDecompLU.h" +#include "TH1D.h" +#include "TFile.h" +#include "TGraph.h" +#include "TCanvas.h" + +//#define DEBUG + +using namespace edm; +using namespace std; + +//---------------------------------------------------------------------------------------------------- + +TH1D *StraightTrackAlignment::newResiduaHist(const char *name) { + return new TH1D(name, ";residual (mm)", 2000, -3., +3.); // in mm +} + +//---------------------------------------------------------------------------------------------------- + +TGraph *newGraph(const string &name, const string &title) { + TGraph *g = new TGraph(); + g->SetName(name.c_str()); + g->SetTitle(title.c_str()); + return g; +} + +//---------------------------------------------------------------------------------------------------- + +StraightTrackAlignment::RPSetPlots::RPSetPlots(const string &_name) : name(_name) { + chisqn_lin_fitted = new TH1D("chi^2 norm, lin, fitted", ";#chi^{2}/ndf;", 5000, 0., 500.); + chisqn_lin_selected = new TH1D("chi^2 norm, lin, selected", ";#chi^{2}/ndf;", 5000, 0., 500.); + chisqn_log_fitted = new TH1D("chi^2 norm, log, fitted", ";log_{10}(#chi^{2}/ndf);", 700, -1., 6.); + chisqn_log_selected = new TH1D("chi^2 norm, log, selected", ";log_{10}(#chi^{2}/ndf);", 700, -1., 6.); + + fitAxVsAyGraph_fitted = newGraph("ax vs. ay, fitted", ";a_{x} (rad);a_{y} (rad)"); + fitAxVsAyGraph_selected = newGraph("ax vs. ay, selected", ";a_{x} (rad);a_{y} (rad)"); + fitBxVsByGraph_fitted = newGraph("bx vs. by, fitted", ";b_{x} (mm);b_{y} (mm)"); + fitBxVsByGraph_selected = newGraph("bx vs. by, selected", ";b_{x} (mm);b_{y} (mm)"); +} +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::RPSetPlots::free() { + delete chisqn_lin_fitted; + delete chisqn_lin_selected; + delete chisqn_log_fitted; + delete chisqn_log_selected; + + delete fitAxVsAyGraph_fitted; + delete fitAxVsAyGraph_selected; + delete fitBxVsByGraph_fitted; + delete fitBxVsByGraph_selected; +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::RPSetPlots::write() const { + chisqn_lin_fitted->Write(); + chisqn_lin_selected->Write(); + chisqn_log_fitted->Write(); + chisqn_log_selected->Write(); + + fitAxVsAyGraph_fitted->Write(); + fitAxVsAyGraph_selected->Write(); + fitBxVsByGraph_fitted->Write(); + fitBxVsByGraph_selected->Write(); +} + +//---------------------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------------------- + +StraightTrackAlignment::StraightTrackAlignment(const ParameterSet &ps) + : verbosity(ps.getUntrackedParameter("verbosity", 0)), + + rpIds(ps.getParameter >("rpIds")), + excludePlanes(ps.getParameter >("excludePlanes")), + z0(ps.getParameter("z0")), + + maxEvents(ps.getParameter("maxEvents")), + + removeImpossible(ps.getParameter("removeImpossible")), + requireNumberOfUnits(ps.getParameter("requireNumberOfUnits")), + requireAtLeast3PotsInOverlap(ps.getParameter("requireAtLeast3PotsInOverlap")), + requireOverlap(ps.getParameter("requireOverlap")), + cutOnChiSqPerNdf(ps.getParameter("cutOnChiSqPerNdf")), + chiSqPerNdfCut(ps.getParameter("chiSqPerNdfCut")), + maxTrackAx(ps.getParameter("maxTrackAx")), + maxTrackAy(ps.getParameter("maxTrackAy")), + + fileNamePrefix(ps.getParameter("fileNamePrefix")), + expandedFileNamePrefix(ps.getParameter("expandedFileNamePrefix")), + factoredFileNamePrefix(ps.getParameter("factoredFileNamePrefix")), + preciseXMLFormat(ps.getParameter("preciseXMLFormat")), + saveXMLUncertainties(ps.getParameter("saveXMLUncertainties")), + + saveIntermediateResults(ps.getParameter("saveIntermediateResults")), + taskDataFileName(ps.getParameter("taskDataFileName")), + taskDataFile(nullptr), + + task(ps), + fitter(ps), + + buildDiagnosticPlots(ps.getParameter("buildDiagnosticPlots")), + diagnosticsFile(ps.getParameter("diagnosticsFile")), + fitNdfHist_fitted(new TH1D("ndf_fitted", ";ndf;", 41, -4.5, 36.5)), + fitNdfHist_selected(new TH1D("ndf_selected", ";ndf;", 41, -4.5, 36.5)), + fitPHist_fitted(new TH1D("p_fitted", ";p value;", 100, 0., 1.)), + fitPHist_selected(new TH1D("p_selected", ";p value;", 100, 0., 1.)), + fitAxHist_fitted(new TH1D("ax_fitted", ";a_{x} (rad);", 10000, -0.1, 0.1)), + fitAxHist_selected(new TH1D("ax_selected", ";a_{x} (rad);", 10000, -0.1, 0.1)), + fitAyHist_fitted(new TH1D("ay_fitted", ";a_{y} (rad);", 10000, -0.1, 0.1)), + fitAyHist_selected(new TH1D("ay_selected", ";a_{y} (rad);", 10000, -0.1, 0.1)), + fitBxHist_fitted(new TH1D("bx_fitted", ";b_{x} (mm);", 500, -30., 30.)), + fitBxHist_selected(new TH1D("bx_selected", ";b_{x} (mm);", 500, -30., 30.)), + fitByHist_fitted(new TH1D("by_fitted", ";b_{y} (mm);", 500, -30., 30.)), + fitByHist_selected(new TH1D("by_selected", ";b_{y} (mm);", 500, -30., 30.)), + + globalPlots("global") { + // open task data file + if (!taskDataFileName.empty()) + taskDataFile = new TFile(taskDataFileName.c_str(), "recreate"); + + // instantiate algorithm objects + // (and save them) + vector alNames(ps.getParameter >("algorithms")); + for (unsigned int i = 0; i < alNames.size(); i++) { + AlignmentAlgorithm *a = nullptr; + + if (alNames[i] == "Ideal") { + IdealResult *ir = new IdealResult(ps, &task); + a = ir; + } else if (alNames[i] == "Jan") { + JanAlignmentAlgorithm *jaa = new JanAlignmentAlgorithm(ps, &task); + a = jaa; + } + + if (a) + algorithms.push_back(a); + else + throw cms::Exception("PPS") << "Unknown alignment algorithm `" << alNames[i] << "'."; + } + + // get constraints type + string ct = ps.getParameter("constraintsType"); + if (ct.compare("fixedDetectors") == 0) + constraintsType = ctFixedDetectors; + else if (ct.compare("standard") == 0) + constraintsType = ctStandard; + else + throw cms::Exception("PPS") << "Unknown constraints type `" << ct << "'."; + + // parse additional accepted RP sets + string aars_str = ps.getParameter("additionalAcceptedRPSets"); + + size_t idx_b = 0, idx_e = string::npos; + while (idx_b != string::npos) { + // get one block - portion between successive ";" + idx_e = aars_str.find(';', idx_b); + size_t len = (idx_e == string::npos) ? string::npos : idx_e - idx_b; + string block = aars_str.substr(idx_b, len); + + // process the block + if (!block.empty()) { + set rpSet; + + // isolate bits (= RP ids) + size_t bi_b = 0, bi_e = string::npos; + while (bi_b != string::npos) { + bi_e = block.find(',', bi_b); + size_t bit_len = (bi_e == string::npos) ? string::npos : bi_e - bi_b; + const string &bit = block.substr(bi_b, bit_len); + + unsigned int rp = atoi(bit.c_str()); + rpSet.insert(rp); + + bi_b = (bi_e == string::npos) ? string::npos : bi_e + 1; + } + + additionalAcceptedRPSets.push_back(rpSet); + } + + // move to next block + idx_b = (idx_e == string::npos) ? string::npos : idx_e + 1; + } +} + +//---------------------------------------------------------------------------------------------------- + +StraightTrackAlignment::~StraightTrackAlignment() { + if (taskDataFile) + delete taskDataFile; + + for (vector::iterator it = algorithms.begin(); it != algorithms.end(); ++it) + delete (*it); + + delete fitNdfHist_fitted; + delete fitNdfHist_selected; + delete fitPHist_fitted; + delete fitPHist_selected; + delete fitAxHist_fitted; + delete fitAyHist_fitted; + delete fitAxHist_selected; + delete fitAyHist_selected; + delete fitBxHist_fitted; + delete fitByHist_fitted; + delete fitBxHist_selected; + delete fitByHist_selected; + + globalPlots.free(); + + for (auto &p : rpSetPlots) + p.second.free(); + + for (auto it = residuaHistograms.begin(); it != residuaHistograms.end(); ++it) { + delete it->second.total_fitted; + delete it->second.total_selected; + delete it->second.selected_vs_chiSq; + for (auto sit = it->second.perRPSet_fitted.begin(); sit != it->second.perRPSet_fitted.end(); ++sit) + delete sit->second; + for (auto sit = it->second.perRPSet_selected.begin(); sit != it->second.perRPSet_selected.end(); ++sit) + delete sit->second; + } +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::begin(edm::ESHandle hRealAlignment, + edm::ESHandle hRealGeometry, + edm::ESHandle hMisalignedGeometry) { + // reset counters + eventsTotal = 0; + eventsFitted = 0; + eventsSelected = 0; + fittedTracksPerRPSet.clear(); + selectedTracksPerRPSet.clear(); + selectedHitsPerPlane.clear(); + + // prepare geometry + task.buildGeometry(rpIds, excludePlanes, hRealGeometry.product(), z0, task.geometry); + + // build matrix index mappings + task.buildIndexMaps(); + + // print geometry info + if (verbosity > 1) + task.geometry.print(); + + // save task (including geometry) and fitter + if (taskDataFile) { + taskDataFile->WriteObject(&task, "task"); + taskDataFile->WriteObject(&fitter, "fitter"); + } + + // initiate the algorithms + for (const auto &a : algorithms) + a->begin(hRealGeometry.product(), hMisalignedGeometry.product()); + + // get initial alignments + initialAlignments = *hRealAlignment; +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::processEvent(const edm::EventID &eventId, + const DetSetVector &uvPatternsStrip, + const DetSetVector &hitsDiamond, + const edm::DetSetVector &hitsPixel, + const DetSetVector &tracksPixel) { + eventsTotal++; + + if (verbosity > 9) + printf("\n---------- StraightTrackAlignment::ProcessEvent (%u:%llu)\n", eventId.run(), eventId.event()); + + // -------------------- STEP 1: get hits from selected RPs + + HitCollection hitSelection; + + // strips + for (const auto &pv : uvPatternsStrip) { + const CTPPSDetId detId(pv.detId()); + const unsigned int rpDecId = detId.arm() * 100 + detId.station() * 10 + detId.rp(); + + // skip if RP not selected + if (find(rpIds.begin(), rpIds.end(), rpDecId) == rpIds.end()) + continue; + + // require exactly 1 U and 1 V pattern, i.e. unique U-V association + unsigned int n_U = 0, n_V = 0; + unsigned int idx_U = 0, idx_V = 0; + for (unsigned int pi = 0; pi < pv.size(); pi++) { + const TotemRPUVPattern &pattern = pv[pi]; + + switch (pattern.projection()) { + case TotemRPUVPattern::projU: + n_U++; + idx_U = pi; + break; + + case TotemRPUVPattern::projV: + n_V++; + idx_V = pi; + break; + + default: + break; + } + } + + if (n_U != 1 || n_V != 1) + continue; + + // skip if patterns not reasonable + if (!pv[idx_U].fittable() || !pv[idx_V].fittable()) + continue; + + // add hits from the U and V pattern + for (const auto &pattern : {pv[idx_U], pv[idx_V]}) { + for (const auto &hitsDetSet : pattern.hits()) { + // skip if sensor not in geometry + if (!task.geometry.isValidSensorId(hitsDetSet.detId())) + continue; + + const double &z = task.geometry.get(hitsDetSet.detId()).z; + for (auto &hit : hitsDetSet) + hitSelection.emplace_back(Hit(hitsDetSet.detId(), 2, hit.position(), hit.sigma(), z)); + } + } + } + + // diamonds + for (const auto &pv : hitsDiamond) { + // skip if RP not selected + const CTPPSDetId detId(pv.detId()); + const unsigned int rpDecId = detId.arm() * 100 + detId.station() * 10 + detId.rp(); + + // skip if RP not selected + if (find(rpIds.begin(), rpIds.end(), rpDecId) == rpIds.end()) + continue; + + // skip if sensor not in geometry + if (!task.geometry.isValidSensorId(detId)) + continue; + + const double &z = task.geometry.get(pv.detId()).z; + + for (const auto &h : pv) + hitSelection.emplace_back(Hit(pv.detId(), 1, h.x(), h.xWidth() / sqrt(12.), z)); + } + + // pixels: data from rec hits + map pixelPlaneMultiplicity; + for (const auto &pv : hitsPixel) + pixelPlaneMultiplicity[pv.detId()] += pv.size(); + + for (const auto &pv : hitsPixel) { + // skip if RP not selected + const CTPPSDetId detId(pv.detId()); + const unsigned int rpDecId = detId.arm() * 100 + detId.station() * 10 + detId.rp(); + + // skip if RP not selected + if (find(rpIds.begin(), rpIds.end(), rpDecId) == rpIds.end()) + continue; + + // skip if sensor not in geometry + if (!task.geometry.isValidSensorId(detId)) + continue; + + // skip if plane multiplicity greater than 1 + if (pixelPlaneMultiplicity[pv.detId()] > 1) + continue; + + for (const auto &h : pv) { + const auto &dg = task.geometry.get(pv.detId()); + const double dz1 = dg.getDirectionData(1).dz; + const double dz2 = dg.getDirectionData(2).dz; + const double z = dg.z + h.point().x() * dz1 + h.point().y() * dz2; + + double x_unc = sqrt(h.error().xx()); + double y_unc = sqrt(h.error().yy()); + + // TODO: Check this + + if (x_unc <= 0.) + x_unc = 10E-3; // mm + if (y_unc <= 0.) + y_unc = 10E-3; // mm + + hitSelection.emplace_back(Hit(pv.detId(), 1, h.point().x(), x_unc, z)); + hitSelection.emplace_back(Hit(pv.detId(), 2, h.point().y(), y_unc, z)); + } + } + + // pixels: data from tracks + for (const auto &pv : tracksPixel) { + const CTPPSDetId rpId(pv.detId()); + const unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); + + // skip if RP not selected + if (find(rpIds.begin(), rpIds.end(), rpDecId) == rpIds.end()) + continue; + + // skip if more than 1 (valid) track in the RP + unsigned int n_valid_tracks = 0; + for (const auto &tr : pv) { + if (tr.isValid()) + n_valid_tracks++; + } + + if (n_valid_tracks > 1) + continue; + + // go through all valid tracks + for (const auto &tr : pv) { + if (!tr.isValid()) + continue; + + // go through all associated rec hits + for (const auto &hv : tr.hits()) { + const CTPPSPixelDetId senId(hv.detId()); + + // skip if sensor not in geometry + if (!task.geometry.isValidSensorId(senId)) + continue; + + for (const auto &h : hv) { + // skip hit if not used for fit + if (!h.isUsedForFit()) + continue; + + const auto &dg = task.geometry.get(senId); + const double dz1 = dg.getDirectionData(1).dz; + const double dz2 = dg.getDirectionData(2).dz; + const double z = dg.z + h.point().x() * dz1 + h.point().y() * dz2; + + double x_unc = sqrt(h.error().xx()); + double y_unc = sqrt(h.error().yy()); + + if (x_unc <= 0.) + x_unc = 10E-3; // mm + if (y_unc <= 0.) + y_unc = 10E-3; // mm + + hitSelection.emplace_back(Hit(senId, 1, h.point().x(), x_unc, z)); + hitSelection.emplace_back(Hit(senId, 2, h.point().y(), y_unc, z)); + } + } + } + } + + if (hitSelection.empty()) + return; + + // -------------------- STEP 2: fit + outlier rejection + + LocalTrackFit trackFit; + if (!fitter.fit(hitSelection, task.geometry, trackFit)) + return; + + set selectedRPs; + for (const auto &hit : hitSelection) { + CTPPSDetId detId(hit.id); + const unsigned int decRPId = detId.arm() * 100 + detId.station() * 10 + detId.rp(); + selectedRPs.insert(decRPId); + } + + eventsFitted++; + fittedTracksPerRPSet[selectedRPs]++; + + // -------------------- STEP 3: quality checks + + bool top = false, bottom = false, horizontal = false; + unordered_set units; + for (const auto &rp : selectedRPs) { + unsigned int rpIdx = rp % 10; + unsigned int stId = rp / 10; + unsigned int unitId = stId * 10; + if (rpIdx > 2) + unitId++; + + if (rpIdx == 0 || rpIdx == 4) + top = true; + if (rpIdx == 1 || rpIdx == 5) + bottom = true; + if (rpIdx == 2 || rpIdx == 3) + horizontal = true; + + units.insert(unitId); + } + + bool overlap = (top && horizontal) || (bottom && horizontal); + + bool rp_set_accepted = true; + + // impossible signature + if (removeImpossible && top && bottom) + rp_set_accepted = false; + + // cleanliness cuts + if (units.size() < requireNumberOfUnits) + rp_set_accepted = false; + + if (requireOverlap && !overlap) + rp_set_accepted = false; + + if (requireAtLeast3PotsInOverlap && overlap && selectedRPs.size() < 3) + rp_set_accepted = false; + + // is it an additional accepted RP set? + if (find(additionalAcceptedRPSets.begin(), additionalAcceptedRPSets.end(), selectedRPs) != + additionalAcceptedRPSets.end()) + rp_set_accepted = true; + + if (verbosity > 5) + printf("* rp set accepted: %u\n", rp_set_accepted); + + bool selected = rp_set_accepted; + + // too bad chisq + if (cutOnChiSqPerNdf && trackFit.chiSqPerNdf() > chiSqPerNdfCut) + selected = false; + + // parallelity cut + if (fabs(trackFit.ax) > maxTrackAx || fabs(trackFit.ay) > maxTrackAy) + selected = false; + + updateDiagnosticHistograms(hitSelection, selectedRPs, trackFit, selected); + + if (verbosity > 5) + printf("* event selected: %u\n", selected); + + if (!selected) + return; + + // update counters + eventsSelected++; + selectedTracksPerRPSet[selectedRPs]++; + + for (const auto &h : hitSelection) + selectedHitsPerPlane[h.id]++; + + // -------------------- STEP 4: FEED ALGORITHMS + + for (auto &a : algorithms) + a->feed(hitSelection, trackFit); + + // -------------------- STEP 5: ENOUGH TRACKS? + + if (eventsSelected == maxEvents) + throw "StraightTrackAlignment: Number of tracks processed reached maximum"; +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::updateDiagnosticHistograms(const HitCollection &selection, + const set &selectedRPs, + const LocalTrackFit &trackFit, + bool trackSelected) { + if (!buildDiagnosticPlots) + return; + + fitNdfHist_fitted->Fill(trackFit.ndf); + fitPHist_fitted->Fill(trackFit.pValue()); + fitAxHist_fitted->Fill(trackFit.ax); + fitAyHist_fitted->Fill(trackFit.ay); + fitBxHist_fitted->Fill(trackFit.bx); + fitByHist_fitted->Fill(trackFit.by); + + globalPlots.chisqn_lin_fitted->Fill(trackFit.chiSqPerNdf()); + globalPlots.chisqn_log_fitted->Fill(log10(trackFit.chiSqPerNdf())); + globalPlots.fitAxVsAyGraph_fitted->SetPoint(globalPlots.fitAxVsAyGraph_fitted->GetN(), trackFit.ax, trackFit.ay); + globalPlots.fitBxVsByGraph_fitted->SetPoint(globalPlots.fitBxVsByGraph_fitted->GetN(), trackFit.bx, trackFit.by); + + if (trackSelected) { + fitNdfHist_selected->Fill(trackFit.ndf); + fitPHist_selected->Fill(trackFit.pValue()); + fitAxHist_selected->Fill(trackFit.ax); + fitAyHist_selected->Fill(trackFit.ay); + fitBxHist_selected->Fill(trackFit.bx); + fitByHist_selected->Fill(trackFit.by); + + globalPlots.chisqn_lin_selected->Fill(trackFit.chiSqPerNdf()); + globalPlots.chisqn_log_selected->Fill(log10(trackFit.chiSqPerNdf())); + globalPlots.fitAxVsAyGraph_selected->SetPoint( + globalPlots.fitAxVsAyGraph_selected->GetN(), trackFit.ax, trackFit.ay); + globalPlots.fitBxVsByGraph_selected->SetPoint( + globalPlots.fitBxVsByGraph_selected->GetN(), trackFit.bx, trackFit.by); + } + + auto it = rpSetPlots.find(selectedRPs); + if (it == rpSetPlots.end()) + it = rpSetPlots.insert({selectedRPs, RPSetPlots(setToString(selectedRPs))}).first; + + it->second.chisqn_lin_fitted->Fill(trackFit.chiSqPerNdf()); + it->second.chisqn_log_fitted->Fill(log10(trackFit.chiSqPerNdf())); + it->second.fitAxVsAyGraph_fitted->SetPoint(it->second.fitAxVsAyGraph_fitted->GetN(), trackFit.ax, trackFit.ay); + it->second.fitBxVsByGraph_fitted->SetPoint(it->second.fitBxVsByGraph_fitted->GetN(), trackFit.bx, trackFit.by); + + if (trackSelected) { + it->second.chisqn_lin_selected->Fill(trackFit.chiSqPerNdf()); + it->second.chisqn_log_selected->Fill(log10(trackFit.chiSqPerNdf())); + it->second.fitAxVsAyGraph_selected->SetPoint(it->second.fitAxVsAyGraph_selected->GetN(), trackFit.ax, trackFit.ay); + it->second.fitBxVsByGraph_selected->SetPoint(it->second.fitBxVsByGraph_selected->GetN(), trackFit.bx, trackFit.by); + } + + for (const auto &hit : selection) { + unsigned int id = hit.id; + + const DetGeometry &geom = task.geometry.get(id); + const auto dirData = geom.getDirectionData(hit.dirIdx); + + double m = hit.position + dirData.s - (hit.z - geom.z) * dirData.dz; + double x = trackFit.ax * hit.z + trackFit.bx; + double y = trackFit.ay * hit.z + trackFit.by; + double f = x * dirData.dx + y * dirData.dy; + double R = m - f; + + auto it = residuaHistograms.find(id); + if (it == residuaHistograms.end()) { + it = residuaHistograms.insert(pair(id, ResiduaHistogramSet())).first; + char buf[30]; + sprintf(buf, "%u: total_fitted", id); + it->second.total_fitted = newResiduaHist(buf); + sprintf(buf, "%u: total_selected", id); + it->second.total_selected = newResiduaHist(buf); + it->second.selected_vs_chiSq = new TGraph(); + sprintf(buf, "%u: selected_vs_chiSq", id); + it->second.selected_vs_chiSq->SetName(buf); + } + + it->second.total_fitted->Fill(R); + + if (trackSelected) { + it->second.total_selected->Fill(R); + it->second.selected_vs_chiSq->SetPoint(it->second.selected_vs_chiSq->GetN(), trackFit.chiSqPerNdf(), R); + } + + auto sit = it->second.perRPSet_fitted.find(selectedRPs); + if (sit == it->second.perRPSet_fitted.end()) { + char buf[10]; + sprintf(buf, "%u: ", id); + string label = buf; + label += setToString(selectedRPs); + sit = + it->second.perRPSet_fitted.insert(pair, TH1D *>(selectedRPs, newResiduaHist(label.c_str()))) + .first; + } + + sit->second->Fill(R); + + if (trackSelected) { + sit = it->second.perRPSet_selected.find(selectedRPs); + if (sit == it->second.perRPSet_selected.end()) { + char buf[10]; + sprintf(buf, "%u: ", id); + string label = buf; + label += setToString(selectedRPs); + sit = it->second.perRPSet_selected + .insert(pair, TH1D *>(selectedRPs, newResiduaHist(label.c_str()))) + .first; + } + + sit->second->Fill(R); + } + } +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::buildConstraints(vector &constraints) { + constraints.clear(); + + switch (constraintsType) { + case ctFixedDetectors: + task.buildFixedDetectorsConstraints(constraints); + return; + + case ctStandard: + task.buildStandardConstraints(constraints); + return; + } +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::finish() { + // print statistics + if (verbosity) { + printf("----------------------------------------------------------------------------------------------------\n"); + printf("\n>> StraightTrackAlignment::Finish\n"); + printf("\tevents total = %i\n", eventsTotal); + printf("\tevents fitted = %i\n", eventsFitted); + printf("\tevents selected = %i\n", eventsSelected); + + printf("\n* events per RP set:\n"); + printf("%30s %10s%10s\n", "set of RPs", "fitted", "selected"); + + for (auto it = fittedTracksPerRPSet.begin(); it != fittedTracksPerRPSet.end(); ++it) { + const string &label = setToString(it->first); + + auto sit = selectedTracksPerRPSet.find(it->first); + unsigned long sv = (sit == selectedTracksPerRPSet.end()) ? 0 : sit->second; + + printf("%30s :%10lu%10lu\n", label.c_str(), it->second, sv); + } + + if (verbosity >= 2) { + printf("\n* hits per plane:\n"); + for (const auto it : selectedHitsPerPlane) { + printf(" "); + printId(it.first); + printf(" : %u\n", it.second); + } + } + } + + // write diagnostics plots + saveDiagnostics(); + + // run analysis + for (auto a : algorithms) + a->analyze(); + + // build constraints + vector constraints; + buildConstraints(constraints); + + // save constraints + if (taskDataFile) + taskDataFile->WriteObject(&constraints, "constraints"); + + if (verbosity) { + printf("\n>> StraightTrackAlignment::Finish > %lu constraints built\n", constraints.size()); + for (unsigned int i = 0; i < constraints.size(); i++) + printf("\t%s\n", constraints[i].name.c_str()); + } + + // solve + vector > results; + for (auto algorithm : algorithms) { + TDirectory *dir = nullptr; + if (taskDataFile && saveIntermediateResults) + dir = taskDataFile->mkdir((algorithm->getName() + "_data").c_str()); + + results.resize(results.size() + 1); + unsigned int rf = algorithm->solve(constraints, results.back(), dir); + + if (rf) + throw cms::Exception("PPS") << "The Solve method of `" << algorithm->getName() + << "' algorithm has failed (return value " << rf << ")."; + } + + // print results + if (verbosity) { + printf("\n>> StraightTrackAlignment::Finish > Print\n"); + + printLineSeparator(results); + printQuantitiesLine(results); + printAlgorithmsLine(results); + + signed int prevRPId = -1; + + for (const auto &dit : task.geometry.getSensorMap()) { + signed int rpId = CTPPSDetId(dit.first).rpId(); + if (rpId != prevRPId) + printLineSeparator(results); + prevRPId = rpId; + + printId(dit.first); + printf(" ║"); + + for (unsigned int q = 0; q < task.quantityClasses.size(); q++) { + for (unsigned int a = 0; a < results.size(); a++) { + const auto it = results[a].find(dit.first); + if (it == results[a].end()) { + if (algorithms[a]->hasErrorEstimate()) + printf("%19s", "----"); + else + printf("%8s", "----"); + + if (a + 1 == results.size()) + printf("║"); + else + printf("│"); + + continue; + } + + const auto &ac = it->second; + double v = 0., e = 0.; + switch (task.quantityClasses[q]) { + case AlignmentTask::qcShR1: + v = ac.getShR1(); + e = ac.getShR1Unc(); + break; + case AlignmentTask::qcShR2: + v = ac.getShR2(); + e = ac.getShR2Unc(); + break; + case AlignmentTask::qcShZ: + v = ac.getShZ(); + e = ac.getShZUnc(); + break; + case AlignmentTask::qcRotZ: + v = ac.getRotZ(); + e = ac.getRotZUnc(); + break; + } + + if (algorithms[a]->hasErrorEstimate()) + printf("%+8.1f ± %7.1f", v * 1E3, e * 1E3); + else + printf("%+8.1f", v * 1E3); + + if (a + 1 == results.size()) + printf("║"); + else + printf("│"); + } + } + + printf("\n"); + } + + printLineSeparator(results); + printAlgorithmsLine(results); + printQuantitiesLine(results); + printLineSeparator(results); + } + + // save results as alignment XML files + for (unsigned int a = 0; a < results.size(); a++) { + // convert readout-direction corrections to X and Y + CTPPSRPAlignmentCorrectionsData convertedAlignments; + for (const auto &p : results[a]) { + const DetGeometry &d = task.geometry.get(p.first); + const auto dir1 = d.getDirectionData(1); + const auto dir2 = d.getDirectionData(2); + + const double det = dir1.dx * dir2.dy - dir1.dy * dir2.dx; + const double sh_x = (dir2.dy * p.second.getShR1() - dir1.dy * p.second.getShR2()) / det; + const double sh_y = (-dir2.dx * p.second.getShR1() + dir1.dx * p.second.getShR2()) / det; + + const double sh_x_e = + sqrt(pow(dir2.dy / det * p.second.getShR1Unc(), 2) + pow(dir1.dy / det * p.second.getShR2Unc(), 2)); + const double sh_y_e = + sqrt(pow(dir2.dx / det * p.second.getShR1Unc(), 2) + pow(dir1.dx / det * p.second.getShR2Unc(), 2)); + + const CTPPSRPAlignmentCorrectionData corr(sh_x, + sh_x_e, + sh_y, + sh_y_e, + p.second.getShZ(), + p.second.getShZUnc(), + 0., + 0., + 0., + 0., + p.second.getRotZ(), + p.second.getRotZUnc()); + convertedAlignments.setSensorCorrection(p.first, corr); + } + + // write non-cumulative alignments + if (!fileNamePrefix.empty()) { + CTPPSRPAlignmentCorrectionsMethods::writeToXML(convertedAlignments, + fileNamePrefix + algorithms[a]->getName() + ".xml", + preciseXMLFormat, + saveXMLUncertainties, + true, + true, + true, + true); + } + + // merge alignments + CTPPSRPAlignmentCorrectionsData cumulativeAlignments; + + cumulativeAlignments.addCorrections(initialAlignments, false, true, true); + cumulativeAlignments.addCorrections(convertedAlignments, false, true, true); + + // write expanded and factored results + if (!expandedFileNamePrefix.empty() || !factoredFileNamePrefix.empty()) { + CTPPSRPAlignmentCorrectionsData expandedAlignments; + CTPPSRPAlignmentCorrectionsData factoredAlignments; + + if (verbosity) + printf(">> Factorizing results of %s algorithm\n", algorithms[a]->getName().c_str()); + + const bool equalWeights = false; + factorRPFromSensorCorrections( + cumulativeAlignments, expandedAlignments, factoredAlignments, task.geometry, equalWeights, verbosity); + + if (!expandedFileNamePrefix.empty()) { + CTPPSRPAlignmentCorrectionsMethods::writeToXML(expandedAlignments, + expandedFileNamePrefix + algorithms[a]->getName() + ".xml", + preciseXMLFormat, + saveXMLUncertainties, + true, + true, + true, + true); + } + + if (!factoredFileNamePrefix.empty()) { + CTPPSRPAlignmentCorrectionsMethods::writeToXML(factoredAlignments, + factoredFileNamePrefix + algorithms[a]->getName() + ".xml", + preciseXMLFormat, + saveXMLUncertainties, + true, + true, + true, + true); + } + } + } + + // prepare algorithms for destructions + for (const auto &algorithm : algorithms) + algorithm->end(); +} + +//---------------------------------------------------------------------------------------------------- + +string StraightTrackAlignment::setToString(const set &s) { + unsigned int N = s.size(); + if (N == 0) + return "empty"; + + string str; + char buf[10]; + unsigned int i = 0; + for (set::iterator it = s.begin(); it != s.end(); ++it, ++i) { + sprintf(buf, "%u", *it); + str += buf; + if (i < N - 1) + str += ", "; + } + + return str; +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::printN(const char *str, unsigned int N) { + for (unsigned int i = 0; i < N; i++) + printf("%s", str); +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::printLineSeparator(const std::vector > &results) { + printf("═════════════════════════╬"); + for (unsigned int q = 0; q < task.quantityClasses.size(); q++) { + for (unsigned int a = 0; a < results.size(); a++) { + printN("═", algorithms[a]->hasErrorEstimate() ? 18 : 8); + if (a + 1 != results.size()) + printf("═"); + } + printf("╬"); + } + printf("\n"); +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::printQuantitiesLine(const std::vector > &results) { + printf(" ║"); + + for (unsigned int q = 0; q < task.quantityClasses.size(); q++) { + unsigned int size = 0; + for (unsigned int a = 0; a < results.size(); a++) + size += (algorithms[a]->hasErrorEstimate()) ? 18 : 8; + size += algorithms.size() - 1; + + const string &tag = task.quantityClassTag(task.quantityClasses[q]); + unsigned int space = (size - tag.size()) / 2; + printN(" ", space); + printf("%s", tag.c_str()); + printN(" ", size - space - tag.size()); + printf("║"); + } + printf("\n"); +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::printAlgorithmsLine(const std::vector > &results) { + printf(" ║"); + + for (unsigned int q = 0; q < task.quantityClasses.size(); q++) { + for (unsigned int a = 0; a < results.size(); a++) { + printf((algorithms[a]->hasErrorEstimate()) ? "%18s" : "%8s", algorithms[a]->getName().substr(0, 8).c_str()); + + if (a + 1 == results.size()) + printf("║"); + else + printf("│"); + } + } + printf("\n"); +} + +//---------------------------------------------------------------------------------------------------- + +void StraightTrackAlignment::saveDiagnostics() const { + if (diagnosticsFile.empty()) + return; + + TFile *df = new TFile(diagnosticsFile.c_str(), "recreate"); + if (df->IsZombie()) + throw cms::Exception("PPS") << "Cannot open file `" << diagnosticsFile << "' for writing."; + + if (buildDiagnosticPlots) { + TDirectory *commonDir = df->mkdir("common"); + gDirectory = commonDir; + + fitNdfHist_fitted->Write(); + fitNdfHist_selected->Write(); + fitAxHist_fitted->Write(); + fitAyHist_fitted->Write(); + fitAxHist_selected->Write(); + fitAyHist_selected->Write(); + fitBxHist_fitted->Write(); + fitByHist_fitted->Write(); + fitBxHist_selected->Write(); + fitByHist_selected->Write(); + fitPHist_fitted->Write(); + fitPHist_selected->Write(); + + gDirectory = commonDir->mkdir("plots global"); + globalPlots.write(); + + TDirectory *ppsDir = commonDir->mkdir("plots per RP set"); + for (map, RPSetPlots>::const_iterator it = rpSetPlots.begin(); it != rpSetPlots.end(); ++it) { + gDirectory = ppsDir->mkdir(setToString(it->first).c_str()); + + it->second.write(); + } + + TDirectory *resDir = commonDir->mkdir("residuals"); + for (map::const_iterator it = residuaHistograms.begin(); + it != residuaHistograms.end(); + ++it) { + char buf[10]; + sprintf(buf, "%u", it->first); + gDirectory = resDir->mkdir(buf); + it->second.total_fitted->Write(); + it->second.total_selected->Write(); + it->second.selected_vs_chiSq->Write(); + + /* + gDirectory = gDirectory->mkdir("fitted per RP set"); + for (map< set, TH1D* >::iterator sit = it->second.perRPSet_fitted.begin(); + sit != it->second.perRPSet_fitted.end(); ++sit) + sit->second->Write(); + gDirectory->cd(".."); +*/ + + gDirectory = gDirectory->mkdir("selected per RP set"); + TCanvas *c = new TCanvas; + c->SetName("alltogether"); + unsigned int idx = 0; + for (map, TH1D *>::const_iterator sit = it->second.perRPSet_selected.begin(); + sit != it->second.perRPSet_selected.end(); + ++sit, ++idx) { + sit->second->SetLineColor(idx + 1); + sit->second->Draw((idx == 0) ? "" : "same"); + sit->second->Write(); + } + c->Write(); + } + } + + // save diagnostics of algorithms + for (vector::const_iterator it = algorithms.begin(); it != algorithms.end(); ++it) { + TDirectory *algDir = df->mkdir((*it)->getName().c_str()); + (*it)->saveDiagnostics(algDir); + } + + delete df; +} diff --git a/CalibPPS/AlignmentRelative/src/Utilities.cc b/CalibPPS/AlignmentRelative/src/Utilities.cc new file mode 100644 index 0000000000000..752a6589427c3 --- /dev/null +++ b/CalibPPS/AlignmentRelative/src/Utilities.cc @@ -0,0 +1,316 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include "CalibPPS/AlignmentRelative/interface/Utilities.h" + +#include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h" +#include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h" +#include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" +#include "DataFormats/CTPPSDetId/interface/CTPPSPixelDetId.h" + +#include "CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h" +#include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsData.h" + +#include "TVectorD.h" +#include "TMatrixD.h" + +#include +#include + +using namespace std; + +//---------------------------------------------------------------------------------------------------- + +void printId(unsigned int id) { + CTPPSDetId detId(id); + + if (detId.subdetId() == CTPPSDetId::sdTrackingStrip) { + TotemRPDetId stDetId(id); + printf("strip %u (%3u.%u)", id, 100 * stDetId.arm() + 10 * stDetId.station() + stDetId.rp(), stDetId.plane()); + } + + if (detId.subdetId() == CTPPSDetId::sdTimingDiamond) { + CTPPSDiamondDetId diDetId(id); + printf("dimnd %u (%3u.%u)", id, 100 * diDetId.arm() + 10 * diDetId.station() + diDetId.rp(), diDetId.plane()); + } + + if (detId.subdetId() == CTPPSDetId::sdTrackingPixel) { + CTPPSPixelDetId piDetId(id); + printf("pixel %u (%3u.%u)", id, 100 * piDetId.arm() + 10 * piDetId.station() + piDetId.rp(), piDetId.plane()); + } +} + +//---------------------------------------------------------------------------------------------------- + +void print(TMatrixD &m, const char *label, bool mathematicaFormat) { + if (mathematicaFormat) { + printf("{"); + for (int i = 0; i < m.GetNrows(); i++) { + if (i > 0) + printf(", "); + printf("{"); + for (int j = 0; j < m.GetNcols(); j++) { + if (j > 0) + printf(", "); + printf("%.3f", m[i][j]); + } + printf("}"); + } + printf("}\n"); + return; + } + + if (label) + printf("\n%s\n", label); + + printf(" | "); + for (int j = 0; j < m.GetNcols(); j++) + printf(" %9i", j); + printf("\n------"); + for (int j = 0; j < m.GetNcols(); j++) + printf("----------"); + printf("\n"); + + for (int i = 0; i < m.GetNrows(); i++) { + printf("%3i | ", i); + for (int j = 0; j < m.GetNcols(); j++) { + double v = m[i][j]; + if (fabs(v) >= 1E4) + printf(" %+9.2E", v); + else if (fabs(v) > 1E-6) + printf(" %+9.2E", v); + else + printf(" 0"); + } + printf("\n"); + } +} + +//---------------------------------------------------------------------------------------------------- + +void factorRPFromSensorCorrections(const CTPPSRPAlignmentCorrectionsData &inputAlignments, + CTPPSRPAlignmentCorrectionsData &expandedAlignments, + CTPPSRPAlignmentCorrectionsData &factoredAlignments, + const AlignmentGeometry &geometry, + bool equalWeights, + unsigned int verbosity) { + // clean first + expandedAlignments.clear(); + factoredAlignments.clear(); + + // save full sensor alignments + map fullAlignments; + map > sensorsPerRP; + for (auto it : inputAlignments.getSensorMap()) { + const auto &sensorId = it.first; + + // with useRPErrors=false the only the sensor uncertainties (coming from the last analysis step) will be used + fullAlignments[sensorId] = inputAlignments.getFullSensorCorrection(sensorId, false); + + sensorsPerRP[CTPPSDetId(sensorId).rpId()].insert(sensorId); + } + + // convert full alignments to expandedAlignments + for (const auto &it : fullAlignments) { + expandedAlignments.setSensorCorrection(it.first, it.second); + } + + // do the factorization RP per RP + for (const auto &rpit : sensorsPerRP) { + CTPPSDetId rpId(rpit.first); + const set &sensors = rpit.second; + + if (verbosity) + printf("* processing RP %u (%u)\n", rpit.first, 100 * rpId.arm() + 10 * rpId.station() + rpId.rp()); + + // determine number of constraints + unsigned int n_constraints = 0; + for (const auto &senId : sensors) { + CTPPSDetId detId(senId); + + if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip) + n_constraints += 1; + + if (rpId.subdetId() == CTPPSDetId::sdTrackingPixel) + n_constraints += 2; + } + + // build matrices + TMatrixD B(n_constraints, 2), Vi(n_constraints, n_constraints), VarM(n_constraints, n_constraints); + TVectorD M(n_constraints); + + double sw2_sh_z = 0., svw2_sh_z = 0., su2w4_sh_z = 0.; + double sw2_rot_x = 0., svw2_rot_x = 0., su2w4_rot_x = 0.; + double sw2_rot_y = 0., svw2_rot_y = 0., su2w4_rot_y = 0.; + double sw2_rot_z = 0., svw2_rot_z = 0., su2w4_rot_z = 0.; + + unsigned int idx = 0; + for (const auto &senId : sensors) { + CTPPSDetId detId(senId); + + const double v_sh_z = fullAlignments[senId].getShZ(); + const double u_sh_z = (fullAlignments[senId].getShZUnc() > 0.) ? fullAlignments[senId].getShZUnc() : 1.; + const double w_sh_z = (equalWeights) ? 1. : 1. / u_sh_z; + sw2_sh_z += w_sh_z * w_sh_z; + svw2_sh_z += v_sh_z * w_sh_z * w_sh_z; + su2w4_sh_z += u_sh_z * u_sh_z * w_sh_z * w_sh_z * w_sh_z * w_sh_z; + + const double v_rot_x = fullAlignments[senId].getRotX(); + const double u_rot_x = (fullAlignments[senId].getRotXUnc() > 0.) ? fullAlignments[senId].getRotXUnc() : 1.; + const double w_rot_x = (equalWeights) ? 1. : 1. / u_rot_x; + sw2_rot_x += w_rot_x * w_rot_x; + svw2_rot_x += v_rot_x * w_rot_x * w_rot_x; + su2w4_rot_x += u_rot_x * u_rot_x * w_rot_x * w_rot_x * w_rot_x * w_rot_x; + + const double v_rot_y = fullAlignments[senId].getRotY(); + const double u_rot_y = (fullAlignments[senId].getRotYUnc() > 0.) ? fullAlignments[senId].getRotYUnc() : 1.; + const double w_rot_y = (equalWeights) ? 1. : 1. / u_rot_y; + sw2_rot_y += w_rot_y * w_rot_y; + svw2_rot_y += v_rot_y * w_rot_y * w_rot_y; + su2w4_rot_y += u_rot_y * u_rot_y * w_rot_y * w_rot_y * w_rot_y * w_rot_y; + + const double v_rot_z = fullAlignments[senId].getRotZ(); + const double u_rot_z = (fullAlignments[senId].getRotZUnc() > 0.) ? fullAlignments[senId].getRotZUnc() : 1.; + const double w_rot_z = (equalWeights) ? 1. : 1. / u_rot_z; + sw2_rot_z += w_rot_z * w_rot_z; + svw2_rot_z += v_rot_z * w_rot_z * w_rot_z; + su2w4_rot_z += u_rot_z * u_rot_z * w_rot_z * w_rot_z * w_rot_z * w_rot_z; + + if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip) { + auto d2 = geometry.get(senId).getDirectionData(2); + + B(idx, 0) = d2.dx; + B(idx, 1) = d2.dy; + + M(idx) = d2.dx * fullAlignments[senId].getShX() + d2.dy * fullAlignments[senId].getShY(); + double unc = + sqrt(pow(d2.dx * fullAlignments[senId].getShXUnc(), 2) + pow(d2.dy * fullAlignments[senId].getShYUnc(), 2)); + if (unc <= 0.) + unc = 1.; + + Vi(idx, idx) = (equalWeights) ? 1. : 1. / unc / unc; + VarM(idx, idx) = unc * unc; + + idx += 1; + } + + if (rpId.subdetId() == CTPPSDetId::sdTrackingPixel) { + B(idx + 0, 0) = 1.; + B(idx + 0, 1) = 0.; + M(idx + 0) = fullAlignments[senId].getShX(); + double x_unc = fullAlignments[senId].getShXUnc(); + if (x_unc <= 0.) + x_unc = 1.; + Vi(idx + 0, idx + 0) = (equalWeights) ? 1. : 1. / x_unc / x_unc; + VarM(idx + 0, idx + 0) = x_unc * x_unc; + + B(idx + 1, 0) = 0.; + B(idx + 1, 1) = 1.; + M(idx + 1) = fullAlignments[senId].getShY(); + double y_unc = fullAlignments[senId].getShYUnc(); + if (y_unc <= 0.) + y_unc = 1.; + Vi(idx + 1, idx + 1) = (equalWeights) ? 1. : 1. / y_unc / y_unc; + VarM(idx + 1, idx + 1) = y_unc * y_unc; + + idx += 2; + } + } + + // calculate per-RP alignment + TMatrixD BT(TMatrixD::kTransposed, B); + TMatrixD BTViB(BT, TMatrixD::kMult, Vi * B); + TMatrixD BTViBi(TMatrixD::kInverted, BTViB); + TMatrixD S(BTViBi * BT * Vi); + TMatrixD ST(TMatrixD::kTransposed, S); + TVectorD th_B(2); + th_B = S * M; + TMatrixD VarTh_B(S * VarM * ST); + + const double m_sh_x = th_B[0], m_sh_x_unc = sqrt(VarTh_B(0, 0)); + const double m_sh_y = th_B[1], m_sh_y_unc = sqrt(VarTh_B(1, 1)); + const double m_sh_z = svw2_sh_z / sw2_sh_z, m_sh_z_unc = sqrt(su2w4_sh_z) / sw2_sh_z; + + const double m_rot_x = svw2_rot_x / sw2_rot_x, m_rot_x_unc = sqrt(su2w4_rot_x) / sw2_rot_x; + const double m_rot_y = svw2_rot_y / sw2_rot_y, m_rot_y_unc = sqrt(su2w4_rot_y) / sw2_rot_y; + const double m_rot_z = svw2_rot_z / sw2_rot_z, m_rot_z_unc = sqrt(su2w4_rot_z) / sw2_rot_z; + + if (verbosity) { + printf(" m_sh_x = (%.1f +- %.1f) um, m_sh_y = (%.1f +- %.1f) um, m_sh_z = (%.1f +- %.1f) mm\n", + m_sh_x * 1E3, + m_sh_x_unc * 1E3, + m_sh_y * 1E3, + m_sh_y_unc * 1E3, + m_sh_z, + m_sh_z_unc); + printf(" m_rot_x = (%.1f +- %.1f) mrad, m_rot_y = (%.1f +- %.1f) mrad, m_rot_z = (%.1f +- %.1f) mrad\n", + m_rot_x * 1E3, + m_rot_x_unc * 1E3, + m_rot_y * 1E3, + m_rot_y_unc * 1E3, + m_rot_z * 1E3, + m_rot_z_unc * 1E3); + } + + factoredAlignments.addRPCorrection(rpId, + CTPPSRPAlignmentCorrectionData(m_sh_x, + m_sh_x_unc, + m_sh_y, + m_sh_y_unc, + m_sh_z, + m_sh_z_unc, + m_rot_x, + m_rot_x_unc, + m_rot_y, + m_rot_y_unc, + m_rot_z, + m_rot_z_unc)); + + // calculate residuals + for (const auto &senId : sensors) { + CTPPSRPAlignmentCorrectionData rc; + + if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip) { + auto d2 = geometry.get(senId).getDirectionData(2); + + const double de_s = + d2.dx * (fullAlignments[senId].getShX() - m_sh_x) + d2.dy * (fullAlignments[senId].getShY() - m_sh_y); + const double de_s_unc = + std::abs(d2.dx * fullAlignments[senId].getShXUnc() + + d2.dy * fullAlignments[senId].getShYUnc()); // the x and y components are fully correlated + + rc = CTPPSRPAlignmentCorrectionData(d2.dx * de_s, + d2.dx * de_s_unc, + d2.dy * de_s, + d2.dy * de_s_unc, + fullAlignments[senId].getShZ() - m_sh_z, + fullAlignments[senId].getShZUnc(), + fullAlignments[senId].getRotX() - m_rot_x, + fullAlignments[senId].getRotXUnc(), + fullAlignments[senId].getRotY() - m_rot_y, + fullAlignments[senId].getRotYUnc(), + fullAlignments[senId].getRotZ() - m_rot_z, + fullAlignments[senId].getRotZUnc()); + } + + if (rpId.subdetId() == CTPPSDetId::sdTrackingPixel) { + rc = CTPPSRPAlignmentCorrectionData(fullAlignments[senId].getShX() - m_sh_x, + fullAlignments[senId].getShXUnc(), + fullAlignments[senId].getShY() - m_sh_y, + fullAlignments[senId].getShYUnc(), + fullAlignments[senId].getShZ() - m_sh_z, + fullAlignments[senId].getShZUnc(), + fullAlignments[senId].getRotX() - m_rot_x, + fullAlignments[senId].getRotXUnc(), + fullAlignments[senId].getRotY() - m_rot_y, + fullAlignments[senId].getRotYUnc(), + fullAlignments[senId].getRotZ() - m_rot_z, + fullAlignments[senId].getRotZUnc()); + } + + factoredAlignments.addSensorCorrection(senId, rc); + } + } +} diff --git a/CalibPPS/AlignmentRelative/src/classes.h b/CalibPPS/AlignmentRelative/src/classes.h new file mode 100644 index 0000000000000..965dfa793e013 --- /dev/null +++ b/CalibPPS/AlignmentRelative/src/classes.h @@ -0,0 +1,10 @@ +#include "CalibPPS/AlignmentRelative/interface/LocalTrackFit.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentTask.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentConstraint.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h" +#include "CalibPPS/AlignmentRelative/interface/JanAlignmentAlgorithm.h" +#include "CalibPPS/AlignmentRelative/interface/AlignmentResult.h" + +#include +#include +#include \ No newline at end of file diff --git a/CalibPPS/AlignmentRelative/src/classes_def.xml b/CalibPPS/AlignmentRelative/src/classes_def.xml new file mode 100644 index 0000000000000..0f5126758a02b --- /dev/null +++ b/CalibPPS/AlignmentRelative/src/classes_def.xml @@ -0,0 +1,36 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/CalibPPS/AlignmentRelative/test/BuildFile.xml b/CalibPPS/AlignmentRelative/test/BuildFile.xml new file mode 100644 index 0000000000000..e464bddcad010 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/BuildFile.xml @@ -0,0 +1,4 @@ + + + + diff --git a/CalibPPS/AlignmentRelative/test/ppsTrackBasedAligAnalRes.cc b/CalibPPS/AlignmentRelative/test/ppsTrackBasedAligAnalRes.cc new file mode 100644 index 0000000000000..fc7a591265d42 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/ppsTrackBasedAligAnalRes.cc @@ -0,0 +1,745 @@ +/**************************************************************************** +* Authors: +* Jan Kašpar (jan.kaspar@gmail.com) +****************************************************************************/ + +#include +#include +#include + +#include "TFile.h" +#include "TGraph.h" +#include "TGraphErrors.h" + +#include "FWCore/Utilities/interface/Exception.h" + +#include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsData.h" +#include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsMethods.h" + +using namespace std; + +//---------------------------------------------------------------------------------------------------- + +TDirectory *mkdir(TDirectory *parent, const char *child) { + TDirectory *dir = (TDirectory *)parent->Get(child); + if (!dir) + dir = parent->mkdir(child); + return dir; +} + +//---------------------------------------------------------------------------------------------------- + +struct Stat { + double s1, sx, sxx, sxxx, sxxxx; + + Stat() : s1(0.), sx(0.), sxx(0.), sxxx(0.), sxxxx(0.) {} + + void fill(double x) { + s1 += 1; + sx += x; + sxx += x * x; + sxxx += x * x * x; + sxxxx += x * x * x * x; + } + + double n() const { return s1; } + + double mean() const { return sx / s1; } + + double meanError() const { + double v = rms() / sqrt(s1); + return v; + } + + double rms() const { + double sig_sq = (sxx - sx * sx / s1) / (s1 - 1.); + if (sig_sq < 0) + sig_sq = 0; + + return sqrt(sig_sq); + } + + double rmsError() const { + // see R.J.Barlow: Statistics, page 78 + double mu = mean(); + double E2 = rms(); + E2 *= E2; + //if (print) printf("\t\t%E, %E, %E, %E, %E\n", sxxxx, -4.*sxxx*mu, 6.*sxx*mu*mu, -4.*sx*mu*mu*mu, s1*mu*mu*mu*mu); + double E4 = (sxxxx - 4. * sxxx * mu + 6. * sxx * mu * mu - 4. * sx * mu * mu * mu + s1 * mu * mu * mu * mu) / s1; + //if (print) printf("\t\tmu = %E, E2 = %E, E4 = %E\n", mu, E2, E4); + double v_sig_sq = (s1 - 1) * ((s1 - 1) * E4 - (s1 - 3) * E2 * E2) / s1 / s1 / s1; + double v_sig = v_sig_sq / 4. / E2; + //if (print) printf("\t\tv_sig_sq = %E\n", v_sig_sq); + if (v_sig < 0 || std::isnan(v_sig) || std::isinf(v_sig)) { + //printf(">> Stat::rmsError > ERROR: v_siq = %E < 0.\n", v_sig); + v_sig = 0; + } + double v = sqrt(v_sig); + return v; + } +}; + +//---------------------------------------------------------------------------------------------------- + +struct DetStat { + // _v ... value given by Jan algorithm + // _u ... uncertainty estimated by Jan algorithm + // _e ... difference (Jan - Ideal) algorithm + Stat sh_x_v, sh_x_u, sh_x_e; + Stat sh_y_v, sh_y_u, sh_y_e; + Stat rot_z_v, rot_z_u, rot_z_e; + Stat sh_z_v, sh_z_u, sh_z_e; +}; + +//---------------------------------------------------------------------------------------------------- + +struct RPStat { + Stat sh_x_v, sh_x_u, sh_x_e; + Stat sh_y_v, sh_y_u, sh_y_e; + Stat rot_z_v, rot_z_u, rot_z_e; + Stat sh_z_v, sh_z_u, sh_z_e; +}; + +//---------------------------------------------------------------------------------------------------- + +struct Desc { + unsigned int N; // number of events + unsigned int i; // iteration + unsigned int d; // detector + Desc(unsigned int _N, unsigned int _i, unsigned int _d) : N(_N), i(_i), d(_d) {} + bool operator<(const Desc &c) const; +}; + +bool Desc::operator<(const Desc &c) const { + if (this->N < c.N) + return true; + if (this->N > c.N) + return false; + if (this->i < c.i) + return true; + if (this->i > c.i) + return false; + if (this->d < c.d) + return true; + if (this->d > c.d) + return false; + return false; +} + +//---------------------------------------------------------------------------------------------------- + +map det_stat; +map rp_stat; + +//---------------------------------------------------------------------------------------------------- + +void resetStatistics() { + det_stat.clear(); + rp_stat.clear(); +} + +//---------------------------------------------------------------------------------------------------- + +void updateSensorStatistics(unsigned int n_events, + unsigned iteration, + const CTPPSRPAlignmentCorrectionsData &r_actual, + const CTPPSRPAlignmentCorrectionsData &r_ideal) { + for (CTPPSRPAlignmentCorrectionsData::mapType::const_iterator it = r_actual.getSensorMap().begin(); + it != r_actual.getSensorMap().end(); + ++it) { + unsigned int id = it->first; + const auto &c_actual = r_actual.getFullSensorCorrection(id, false); + const auto &c_ideal = r_ideal.getFullSensorCorrection(id, false); + + DetStat &s = det_stat[Desc(n_events, iteration, id)]; + + s.sh_x_v.fill(c_actual.getShX()); + s.sh_x_u.fill(c_actual.getShXUnc()); + s.sh_x_e.fill(c_actual.getShX() - c_ideal.getShX()); + + s.sh_y_v.fill(c_actual.getShY()); + s.sh_y_u.fill(c_actual.getShYUnc()); + s.sh_y_e.fill(c_actual.getShY() - c_ideal.getShY()); + + s.sh_z_v.fill(c_actual.getShZ()); + s.sh_z_u.fill(c_actual.getShZUnc()); + s.sh_z_e.fill(c_actual.getShZ() - c_ideal.getShZ()); + + s.rot_z_v.fill(c_actual.getRotZ()); + s.rot_z_u.fill(c_actual.getRotZUnc()); + s.rot_z_e.fill(c_actual.getRotZ() - c_ideal.getRotZ()); + } +} + +//---------------------------------------------------------------------------------------------------- + +void updateRPStatistics(unsigned int n_events, + unsigned iteration, + const CTPPSRPAlignmentCorrectionsData &r_actual, + const CTPPSRPAlignmentCorrectionsData &r_ideal) { + for (CTPPSRPAlignmentCorrectionsData::mapType::const_iterator it = r_actual.getRPMap().begin(); + it != r_actual.getRPMap().end(); + ++it) { + unsigned int id = it->first; + const auto &c_actual = r_actual.getRPCorrection(id); + const auto &c_ideal = r_ideal.getRPCorrection(id); + + RPStat &s = rp_stat[Desc(n_events, iteration, id)]; + + s.sh_x_v.fill(c_actual.getShX()); + s.sh_x_u.fill(c_actual.getShXUnc()); + s.sh_x_e.fill(c_actual.getShX() - c_ideal.getShX()); + + s.sh_y_v.fill(c_actual.getShY()); + s.sh_y_u.fill(c_actual.getShYUnc()); + s.sh_y_e.fill(c_actual.getShY() - c_ideal.getShY()); + + s.sh_z_v.fill(c_actual.getShZ()); + s.sh_z_u.fill(c_actual.getShZUnc()); + s.sh_z_e.fill(c_actual.getShZ() - c_ideal.getShZ()); + + s.rot_z_v.fill(c_actual.getRotZ()); + s.rot_z_u.fill(c_actual.getRotZUnc()); + s.rot_z_e.fill(c_actual.getRotZ() - c_ideal.getRotZ()); + } +} + +//---------------------------------------------------------------------------------------------------- + +struct StatGraphs { + TGraph *n; + TGraphErrors *v_m, *v_v, *u_m, *u_v, *e_m, *e_v, *eR; + StatGraphs() + : n(new TGraph()), + v_m(new TGraphErrors()), + v_v(new TGraphErrors()), + u_m(new TGraphErrors()), + u_v(new TGraphErrors()), + e_m(new TGraphErrors()), + e_v(new TGraphErrors()), + eR(new TGraphErrors()) {} + + void write(const char *xLabel); +}; + +//---------------------------------------------------------------------------------------------------- + +#define ENTRY(tag, label) \ + tag->SetName(#tag); \ + sprintf(buf, ";%s;" label, xLabel); \ + tag->SetTitle(buf); \ + tag->Write(); + +void StatGraphs::write(const char *xLabel) { + char buf[50]; + ENTRY(n, "number of repetitions"); + ENTRY(v_m, "value mean"); + ENTRY(v_v, "value variation"); + ENTRY(u_m, "estim. uncertainty mean"); + ENTRY(u_v, "estim. uncertainty variation"); + ENTRY(e_m, "error mean"); + ENTRY(e_v, "error variation"); + ENTRY(eR, "error variation / estim. uncertainty"); +} + +//---------------------------------------------------------------------------------------------------- + +struct DetGraphs { + StatGraphs sh_x, sh_y, rot_z, sh_z; + void fill(double x, const DetStat &); + void write(const char *xLabel); +}; + +//---------------------------------------------------------------------------------------------------- + +double eR_error(const Stat &e, const Stat &u) { + double a = e.rms(), ae = e.rmsError(); + double b = u.mean(), be = u.meanError(); + + return (b <= 0) ? 0. : a / b * sqrt(ae * ae / a / a + be * be / b / b); +} + +//---------------------------------------------------------------------------------------------------- + +void DetGraphs::fill(double x, const DetStat &s) { + int idx = sh_x.n->GetN(); + + sh_x.n->SetPoint(idx, x, s.sh_x_u.n()); + sh_x.v_m->SetPoint(idx, x, s.sh_x_v.mean()); + sh_x.v_v->SetPoint(idx, x, s.sh_x_v.rms()); + sh_x.u_m->SetPoint(idx, x, s.sh_x_u.mean()); + sh_x.u_v->SetPoint(idx, x, s.sh_x_u.rms()); + sh_x.e_m->SetPoint(idx, x, s.sh_x_e.mean()); + sh_x.e_v->SetPoint(idx, x, s.sh_x_e.rms()); + sh_x.eR->SetPoint(idx, x, s.sh_x_e.rms() / s.sh_x_u.mean()); + + sh_x.v_m->SetPointError(idx, 0., s.sh_x_v.meanError()); + sh_x.v_v->SetPointError(idx, 0., s.sh_x_v.rmsError()); + sh_x.u_m->SetPointError(idx, 0., s.sh_x_u.meanError()); + sh_x.u_v->SetPointError(idx, 0., s.sh_x_u.rmsError()); + sh_x.e_m->SetPointError(idx, 0., s.sh_x_e.meanError()); + sh_x.e_v->SetPointError(idx, 0., s.sh_x_e.rmsError()); + sh_x.eR->SetPointError(idx, 0., eR_error(s.sh_x_e, s.sh_x_u)); + + sh_y.n->SetPoint(idx, x, s.sh_y_u.n()); + sh_y.v_m->SetPoint(idx, x, s.sh_y_v.mean()); + sh_y.v_v->SetPoint(idx, x, s.sh_y_v.rms()); + sh_y.u_m->SetPoint(idx, x, s.sh_y_u.mean()); + sh_y.u_v->SetPoint(idx, x, s.sh_y_u.rms()); + sh_y.e_m->SetPoint(idx, x, s.sh_y_e.mean()); + sh_y.e_v->SetPoint(idx, x, s.sh_y_e.rms()); + sh_y.eR->SetPoint(idx, x, s.sh_y_e.rms() / s.sh_y_u.mean()); + + sh_y.v_m->SetPointError(idx, 0., s.sh_y_v.meanError()); + sh_y.v_v->SetPointError(idx, 0., s.sh_y_v.rmsError()); + sh_y.u_m->SetPointError(idx, 0., s.sh_y_u.meanError()); + sh_y.u_v->SetPointError(idx, 0., s.sh_y_u.rmsError()); + sh_y.e_m->SetPointError(idx, 0., s.sh_y_e.meanError()); + sh_y.e_v->SetPointError(idx, 0., s.sh_y_e.rmsError()); + sh_y.eR->SetPointError(idx, 0., eR_error(s.sh_y_e, s.sh_y_u)); + + rot_z.n->SetPoint(idx, x, s.rot_z_u.n()); + rot_z.v_m->SetPoint(idx, x, s.rot_z_v.mean()); + rot_z.v_v->SetPoint(idx, x, s.rot_z_v.rms()); + rot_z.u_m->SetPoint(idx, x, s.rot_z_u.mean()); + rot_z.u_v->SetPoint(idx, x, s.rot_z_u.rms()); + rot_z.e_m->SetPoint(idx, x, s.rot_z_e.mean()); + rot_z.e_v->SetPoint(idx, x, s.rot_z_e.rms()); + rot_z.eR->SetPoint(idx, x, s.rot_z_e.rms() / s.rot_z_u.mean()); + + rot_z.v_m->SetPointError(idx, 0., s.rot_z_v.meanError()); + rot_z.v_v->SetPointError(idx, 0., s.rot_z_v.rmsError()); + rot_z.u_m->SetPointError(idx, 0., s.rot_z_u.meanError()); + rot_z.u_v->SetPointError(idx, 0., s.rot_z_u.rmsError()); + rot_z.e_m->SetPointError(idx, 0., s.rot_z_e.meanError()); + rot_z.e_v->SetPointError(idx, 0., s.rot_z_e.rmsError()); + rot_z.eR->SetPointError(idx, 0., eR_error(s.rot_z_e, s.rot_z_u)); + + sh_z.n->SetPoint(idx, x, s.sh_z_u.n()); + sh_z.v_m->SetPoint(idx, x, s.sh_z_v.mean()); + sh_z.v_v->SetPoint(idx, x, s.sh_z_v.rms()); + sh_z.u_m->SetPoint(idx, x, s.sh_z_u.mean()); + sh_z.u_v->SetPoint(idx, x, s.sh_z_u.rms()); + sh_z.e_m->SetPoint(idx, x, s.sh_z_e.mean()); + sh_z.e_v->SetPoint(idx, x, s.sh_z_e.rms()); + sh_z.eR->SetPoint(idx, x, s.sh_z_e.rms() / s.sh_z_u.mean()); + + sh_z.v_m->SetPointError(idx, 0., s.sh_z_v.meanError()); + sh_z.v_v->SetPointError(idx, 0., s.sh_z_v.rmsError()); + sh_z.u_m->SetPointError(idx, 0., s.sh_z_u.meanError()); + sh_z.u_v->SetPointError(idx, 0., s.sh_z_u.rmsError()); + sh_z.e_m->SetPointError(idx, 0., s.sh_z_e.meanError()); + sh_z.e_v->SetPointError(idx, 0., s.sh_z_e.rmsError()); + sh_z.eR->SetPointError(idx, 0., eR_error(s.sh_z_e, s.sh_z_u)); +} + +//---------------------------------------------------------------------------------------------------- + +void DetGraphs::write(const char *xLabel) { + gDirectory = mkdir(gDirectory, "sh_x"); + sh_x.write(xLabel); + gDirectory->cd(".."); + gDirectory = mkdir(gDirectory, "sh_y"); + sh_y.write(xLabel); + gDirectory->cd(".."); + gDirectory = mkdir(gDirectory, "rot_z"); + rot_z.write(xLabel); + gDirectory->cd(".."); + gDirectory = mkdir(gDirectory, "sh_z"); + sh_z.write(xLabel); + gDirectory->cd(".."); +} +//---------------------------------------------------------------------------------------------------- + +struct RPGraphs { + StatGraphs sh_x, sh_y, rot_z, sh_z; + void fill(double x, const RPStat &); + void write(const char *xLabel); +}; + +//---------------------------------------------------------------------------------------------------- + +void RPGraphs::fill(double x, const RPStat &s) { + int idx = sh_x.n->GetN(); + + sh_x.n->SetPoint(idx, x, s.sh_x_u.n()); + sh_x.v_m->SetPoint(idx, x, s.sh_x_v.mean()); + sh_x.v_v->SetPoint(idx, x, s.sh_x_v.rms()); + sh_x.u_m->SetPoint(idx, x, s.sh_x_u.mean()); + sh_x.u_v->SetPoint(idx, x, s.sh_x_u.rms()); + sh_x.e_m->SetPoint(idx, x, s.sh_x_e.mean()); + sh_x.e_v->SetPoint(idx, x, s.sh_x_e.rms()); + sh_x.eR->SetPoint(idx, x, s.sh_x_e.rms() / s.sh_x_u.mean()); + + sh_x.v_m->SetPointError(idx, 0., s.sh_x_v.meanError()); + sh_x.v_v->SetPointError(idx, 0., s.sh_x_v.rmsError()); + sh_x.u_m->SetPointError(idx, 0., s.sh_x_u.meanError()); + sh_x.u_v->SetPointError(idx, 0., s.sh_x_u.rmsError()); + sh_x.e_m->SetPointError(idx, 0., s.sh_x_e.meanError()); + sh_x.e_v->SetPointError(idx, 0., s.sh_x_e.rmsError()); + sh_x.eR->SetPointError(idx, 0., eR_error(s.sh_x_e, s.sh_x_u)); + + sh_y.n->SetPoint(idx, x, s.sh_y_u.n()); + sh_y.v_m->SetPoint(idx, x, s.sh_y_v.mean()); + sh_y.v_v->SetPoint(idx, x, s.sh_y_v.rms()); + sh_y.u_m->SetPoint(idx, x, s.sh_y_u.mean()); + sh_y.u_v->SetPoint(idx, x, s.sh_y_u.rms()); + sh_y.e_m->SetPoint(idx, x, s.sh_y_e.mean()); + sh_y.e_v->SetPoint(idx, x, s.sh_y_e.rms()); + sh_y.eR->SetPoint(idx, x, s.sh_y_e.rms() / s.sh_y_u.mean()); + + sh_y.v_m->SetPointError(idx, 0., s.sh_y_v.meanError()); + sh_y.v_v->SetPointError(idx, 0., s.sh_y_v.rmsError()); + sh_y.u_m->SetPointError(idx, 0., s.sh_y_u.meanError()); + sh_y.u_v->SetPointError(idx, 0., s.sh_y_u.rmsError()); + sh_y.e_m->SetPointError(idx, 0., s.sh_y_e.meanError()); + sh_y.e_v->SetPointError(idx, 0., s.sh_y_e.rmsError()); + sh_y.eR->SetPointError(idx, 0., eR_error(s.sh_y_e, s.sh_y_u)); + + rot_z.n->SetPoint(idx, x, s.rot_z_u.n()); + rot_z.v_m->SetPoint(idx, x, s.rot_z_v.mean()); + rot_z.v_v->SetPoint(idx, x, s.rot_z_v.rms()); + rot_z.u_m->SetPoint(idx, x, s.rot_z_u.mean()); + rot_z.u_v->SetPoint(idx, x, s.rot_z_u.rms()); + rot_z.e_m->SetPoint(idx, x, s.rot_z_e.mean()); + rot_z.e_v->SetPoint(idx, x, s.rot_z_e.rms()); + rot_z.eR->SetPoint(idx, x, s.rot_z_e.rms() / s.rot_z_u.mean()); + + rot_z.v_m->SetPointError(idx, 0., s.rot_z_v.meanError()); + rot_z.v_v->SetPointError(idx, 0., s.rot_z_v.rmsError()); + rot_z.u_m->SetPointError(idx, 0., s.rot_z_u.meanError()); + rot_z.u_v->SetPointError(idx, 0., s.rot_z_u.rmsError()); + rot_z.e_m->SetPointError(idx, 0., s.rot_z_e.meanError()); + rot_z.e_v->SetPointError(idx, 0., s.rot_z_e.rmsError()); + rot_z.eR->SetPointError(idx, 0., eR_error(s.rot_z_e, s.rot_z_u)); + + sh_z.n->SetPoint(idx, x, s.sh_z_u.n()); + sh_z.v_m->SetPoint(idx, x, s.sh_z_v.mean()); + sh_z.v_v->SetPoint(idx, x, s.sh_z_v.rms()); + sh_z.u_m->SetPoint(idx, x, s.sh_z_u.mean()); + sh_z.u_v->SetPoint(idx, x, s.sh_z_u.rms()); + sh_z.e_m->SetPoint(idx, x, s.sh_z_e.mean()); + sh_z.e_v->SetPoint(idx, x, s.sh_z_e.rms()); + sh_z.eR->SetPoint(idx, x, s.sh_z_e.rms() / s.sh_z_u.mean()); + + sh_z.v_m->SetPointError(idx, 0., s.sh_z_v.meanError()); + sh_z.v_v->SetPointError(idx, 0., s.sh_z_v.rmsError()); + sh_z.u_m->SetPointError(idx, 0., s.sh_z_u.meanError()); + sh_z.u_v->SetPointError(idx, 0., s.sh_z_u.rmsError()); + sh_z.e_m->SetPointError(idx, 0., s.sh_z_e.meanError()); + sh_z.e_v->SetPointError(idx, 0., s.sh_z_e.rmsError()); + sh_z.eR->SetPointError(idx, 0., eR_error(s.sh_z_e, s.sh_z_u)); +} + +//---------------------------------------------------------------------------------------------------- + +void RPGraphs::write(const char *xLabel) { + gDirectory = mkdir(gDirectory, "sh_x"); + sh_x.write(xLabel); + gDirectory->cd(".."); + gDirectory = mkdir(gDirectory, "sh_y"); + sh_y.write(xLabel); + gDirectory->cd(".."); + gDirectory = mkdir(gDirectory, "rot_z"); + rot_z.write(xLabel); + gDirectory->cd(".."); + gDirectory = mkdir(gDirectory, "sh_z"); + sh_z.write(xLabel); + gDirectory->cd(".."); +} + +//---------------------------------------------------------------------------------------------------- + +TFile *sf; // output file + +void writeGraphs(string r, string o, string d) { + map dFcnN; // here Desc::N is 0 by definition + map dFcnIt; // here Desc::i is 0 by definition + for (map::iterator it = det_stat.begin(); it != det_stat.end(); ++it) { + dFcnN[Desc(0, it->first.i, it->first.d)].fill(it->first.N, it->second); + dFcnIt[Desc(it->first.N, 0, it->first.d)].fill(it->first.i, it->second); + } + + map rFcnN; // here Desc::N is 0 by definition + map rFcnIt; // here Desc::i is 0 by definition + for (map::iterator it = rp_stat.begin(); it != rp_stat.end(); ++it) { + rFcnN[Desc(0, it->first.i, it->first.d)].fill(it->first.N, it->second); + rFcnIt[Desc(it->first.N, 0, it->first.d)].fill(it->first.i, it->second); + } + + r = r.replace(r.find(':'), 1, ">"); + o = o.replace(o.find(':'), 1, ">"); + d = d.replace(d.find(':'), 1, ">"); + + gDirectory = sf; + gDirectory = mkdir(gDirectory, r.c_str()); + gDirectory = mkdir(gDirectory, o.c_str()); + gDirectory = mkdir(gDirectory, d.c_str()); + + char buf[100]; + gDirectory = mkdir(gDirectory, "fcn_of_N"); + for (map::iterator it = dFcnN.begin(); it != dFcnN.end(); ++it) { + sprintf(buf, "iteration>%u", it->first.i); + gDirectory = mkdir(gDirectory, buf); + sprintf(buf, "%u", it->first.d); + gDirectory = mkdir(gDirectory, buf); + it->second.write("tracks"); + gDirectory->cd("../.."); + } + + for (map::iterator it = rFcnN.begin(); it != rFcnN.end(); ++it) { + sprintf(buf, "iteration>%u", it->first.i); + gDirectory = mkdir(gDirectory, buf); + sprintf(buf, "RP %u", it->first.d); + gDirectory = mkdir(gDirectory, buf); + it->second.write("tracks"); + gDirectory->cd("../.."); + } + gDirectory->cd(".."); + + gDirectory = mkdir(gDirectory, "fcn_of_iteration"); + for (map::iterator it = dFcnIt.begin(); it != dFcnIt.end(); ++it) { + sprintf(buf, "N>%u", it->first.N); + gDirectory = mkdir(gDirectory, buf); + sprintf(buf, "%u", it->first.d); + gDirectory = mkdir(gDirectory, buf); + it->second.write("iteration"); + gDirectory->cd("../.."); + } + + for (map::iterator it = rFcnIt.begin(); it != rFcnIt.end(); ++it) { + sprintf(buf, "N>%u", it->first.N); + gDirectory = mkdir(gDirectory, buf); + sprintf(buf, "RP %u", it->first.d); + gDirectory = mkdir(gDirectory, buf); + it->second.write("iteration"); + gDirectory->cd("../.."); + } + gDirectory->cd(".."); +} + +//---------------------------------------------------------------------------------------------------- + +bool isRegDir(const dirent *de) { + if (de->d_type != DT_DIR) + return false; + + if (!strcmp(de->d_name, ".")) + return false; + + if (!strcmp(de->d_name, "..")) + return false; + + return true; +} + +//---------------------------------------------------------------------------------------------------- + +// Choices for RP input. +// * First of all, the Jan to Ideal comparison is not possible for RP data - +// the factorizations are different, because of zero errors of the Ideal +// algorithm. +// * PLAIN results do not contain RP data, do not use. +// * CUMULATIVE contain the factorization from the previous iteration. +// * FACTORED is probably the best choice. +string r_actual_file("./cumulative_factored_results_Jan.xml"); +string r_ideal_file("./cumulative_factored_results_Ideal.xml"); + +// Choices for sensor input. +// * PLAIN results contain only the correction by the last iteration. +// * Expansion is not done, hence CUMULATIVE results contain internal ALIGNMENT +// only. Factorization inherited from the previous iteration. +// * FACTORED results shall not be used - factorization is different for Jan and +// Ideal algorithms (Ideal has zero errors). +// * If EXPANDED results are available, they might be used, they contain both +// internal and RP alignment. Hence comparison between different misalignments +// is not possible. +string s_actual_file("./cumulative_expanded_results_Jan.xml"); +string s_ideal_file("./cumulative_expanded_results_Ideal.xml"); + +//---------------------------------------------------------------------------------------------------- + +CTPPSRPAlignmentCorrectionsData LoadAlignment(const string &fn) { + const auto &seq = CTPPSRPAlignmentCorrectionsMethods::loadFromXML(fn); + if (seq.empty()) + throw cms::Exception("PPS") << "LoadAlignment: alignment sequence empty, file: " << fn; + + return seq[0].second; +} + +//---------------------------------------------------------------------------------------------------- + +void PrintHelp(const char *name) { + printf("USAGE: %s r_actual_file r_ideal_file s_actual_file s_ideal_file\n", name); + printf(" %s --help (to print this help)\n", name); + printf("PARAMETERS:\n"); + printf("\tr_actual_file\t file with actual RP results (default %s)\n", r_actual_file.c_str()); + printf("\tr_ideal_file\t file with ideal RP results (default %s)\n", r_ideal_file.c_str()); + printf("\ts_actual_file\t file with actual sensor results (default %s)\n", s_actual_file.c_str()); + printf("\ts_ideal_file\t file with ideal sensor results (default %s)\n", s_ideal_file.c_str()); +} + +//---------------------------------------------------------------------------------------------------- + +int main(int argc, const char *argv[]) { + if (argc > 1) { + if (!strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) { + PrintHelp(argv[0]); + return 0; + } + + if (argc > 1) + r_actual_file = argv[1]; + if (argc > 2) + r_ideal_file = argv[2]; + if (argc > 3) + s_actual_file = argv[3]; + if (argc > 4) + s_ideal_file = argv[4]; + } + + printf("r_actual_file: %s\n", r_actual_file.c_str()); + printf("r_ideal_file: %s\n", r_ideal_file.c_str()); + printf("s_actual_file: %s\n", s_actual_file.c_str()); + printf("s_ideal_file: %s\n", s_ideal_file.c_str()); + + // open output file + sf = new TFile("result_summary.root", "recreate"); + + // traverse directory structure + try { + DIR *dp_r = opendir("."); + dirent *de_r; + // traverse RPs directories + while ((de_r = readdir(dp_r))) { + if (!isRegDir(de_r)) + continue; + + chdir(de_r->d_name); + DIR *dp_o = opendir("."); + // traverse optimized directories + dirent *de_o; + while ((de_o = readdir(dp_o))) { + if (!isRegDir(de_o)) + continue; + + chdir(de_o->d_name); + DIR *dp_d = opendir("."); + dirent *de_d; + // traverse tr_dist directories + while ((de_d = readdir(dp_d))) { + if (!isRegDir(de_d)) + continue; + + // sort the tr_N directories by N + chdir(de_d->d_name); + DIR *dp_N = opendir("."); + dirent *de_N; + map nMap; + while ((de_N = readdir(dp_N))) { + if (isRegDir(de_N)) { + string sN = de_N->d_name; + sN = sN.substr(sN.find(':') + 1); + unsigned int N = atof(sN.c_str()); + nMap[N] = de_N->d_name; + } + } + + resetStatistics(); + + // traverse tr_N directories + for (map::iterator nit = nMap.begin(); nit != nMap.end(); ++nit) { + chdir(nit->second.c_str()); + DIR *dp_i = opendir("."); + dirent *de_i; + // traverse repetitions directories + while ((de_i = readdir(dp_i))) { + if (!isRegDir(de_i)) + continue; + + chdir(de_i->d_name); + + // sort the iteration directories + DIR *dp_it = opendir("."); + dirent *de_it; + map itMap; + while ((de_it = readdir(dp_it))) { + if (isRegDir(de_it)) { + unsigned idx = 9; + if (de_it->d_name[9] == ':') + idx = 10; + unsigned int it = atoi(de_it->d_name + idx); + itMap[it] = de_it->d_name; + } + } + + // traverse iteration directories + for (map::iterator iit = itMap.begin(); iit != itMap.end(); ++iit) { + chdir(iit->second.c_str()); + + printf("%s|%s|%s|%s|%s|%s\n", + de_r->d_name, + de_o->d_name, + de_d->d_name, + nit->second.c_str(), + de_i->d_name, + iit->second.c_str()); + + // load and process alignments + try { + const auto &r_actual = LoadAlignment(r_actual_file); + const auto &r_ideal = LoadAlignment(r_ideal_file); + updateRPStatistics(nit->first, iit->first, r_actual, r_ideal); + + const auto &s_actual = LoadAlignment(s_actual_file); + const auto &s_ideal = LoadAlignment(s_ideal_file); + updateSensorStatistics(nit->first, iit->first, s_actual, s_ideal); + } catch (cms::Exception &e) { + printf("ERROR: A CMS exception has been caught:\n%s\nSkipping this directory.\n", e.what()); + } + + chdir(".."); + } + + closedir(dp_it); + chdir(".."); + } + + closedir(dp_i); + chdir(".."); + } + + writeGraphs(de_r->d_name, de_o->d_name, de_d->d_name); + + closedir(dp_N); + chdir(".."); + } + + closedir(dp_d); + chdir(".."); + } + + closedir(dp_o); + chdir(".."); + } + + closedir(dp_r); + chdir(".."); + + delete sf; + } + + catch (cms::Exception &e) { + printf("ERROR: A CMS exception has been caught:\n%s\nStopping.\n", e.what()); + } + + catch (std::exception &e) { + printf("ERROR: A std::exception has been caught:\n%s\nStopping.\n", e.what()); + } + + catch (...) { + printf("ERROR: An exception has been caught, stopping.\n"); + } + + return 0; +} diff --git a/CalibPPS/AlignmentRelative/test/run_all_tests b/CalibPPS/AlignmentRelative/test/run_all_tests new file mode 100755 index 0000000000000..396092e1e82a9 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/run_all_tests @@ -0,0 +1,17 @@ +#!/bin/bash + +function RunOneTest() +{ + local dir="$1" + + cd "$dir" + ./run_test +} + +RunOneTest "test_modify_singular_modes" & + +RunOneTest "test_with_data" & + +RunOneTest "test_with_mc/iterations" & +RunOneTest "test_with_mc/simple" & +RunOneTest "test_with_mc/statistics" & \ No newline at end of file diff --git a/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/input.xml b/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/input.xml new file mode 100644 index 0000000000000..56650cfcbd05b --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/input.xml @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/run_test b/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/run_test new file mode 100755 index 0000000000000..e57fbb20dd81e --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/run_test @@ -0,0 +1,24 @@ +#!/bin/bash + +mkdir -p "data" +cd "data" + +tag="test_input" +cat "../template_tba_cfg.py" | sed "\ + s|\$misalignment_file|input.xml|;\ + s|\$output_iteration|${tag}_results_iteration_|;\ + s|\$output_cumulative|${tag}_results_cumulative_|;\ + " > "${tag}_cfg.py" +cmsRun "${tag}_cfg.py" &> "${tag}.log" + +tag="msm" +cp "../template_msm_cfg.py" "${tag}_cfg.py" +cmsRun "${tag}_cfg.py" &> "${tag}.log" + +tag="test_output" +cat "../template_tba_cfg.py" | sed "\ + s|\$misalignment_file|data/output.xml|;\ + s|\$output_iteration|${tag}_results_iteration_|;\ + s|\$output_cumulative|${tag}_results_cumulative_|;\ + " > "${tag}_cfg.py" +cmsRun "${tag}_cfg.py" &> "${tag}.log" diff --git a/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/template_msm_cfg.py b/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/template_msm_cfg.py new file mode 100644 index 0000000000000..665e476b380ff --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/template_msm_cfg.py @@ -0,0 +1,56 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("ppsModifySingularModesTests") + +# minimum of logs +process.MessageLogger = cms.Service("MessageLogger", + statistics = cms.untracked.vstring(), + destinations = cms.untracked.vstring('cout'), + cout = cms.untracked.PSet( + threshold = cms.untracked.string('WARNING') + ) +) + +# data source +process.source = cms.Source("EmptySource") + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(1) +) + +# geometry +process.load("Geometry.VeryForwardGeometry.geometryRPFromDD_2018_cfi") +del(process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles[-1]) +process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles.append("CalibPPS/AlignmentRelative/test/test_with_mc/RP_Dist_Beam_Cent.xml") + +# input +input = "CalibPPS/AlignmentRelative/test/test_modify_singular_modes/input.xml" + +# initial alignments +process.load("CalibPPS.ESProducers.ctppsRPAlignmentCorrectionsDataESSourceXML_cfi") +process.ctppsRPAlignmentCorrectionsDataESSourceXML.RealFiles = cms.vstring(input) + +# worker +process.load("CalibPPS.AlignmentRelative.ppsModifySingularModes_cfi") +process.ppsModifySingularModes.inputFile = input +process.ppsModifySingularModes.outputFile = 'output.xml' + +# x, y and z: mm +# rho in rad + +process.ppsModifySingularModes.z1 = 213.000 * 1E3 # 56-210-fr-vr +process.ppsModifySingularModes.z2 = 220.000 * 1E3 # 56-220-fr-vr + +process.ppsModifySingularModes.de_x1 = +0.100 +process.ppsModifySingularModes.de_x2 = -0.100 + +process.ppsModifySingularModes.de_y1 = +0.200 +process.ppsModifySingularModes.de_y2 = -0.200 + +process.ppsModifySingularModes.de_rho1 = -0.005 +process.ppsModifySingularModes.de_rho2 = -0.010 + +# processing sequence +process.p = cms.Path( + process.ppsModifySingularModes +) diff --git a/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/template_tba_cfg.py b/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/template_tba_cfg.py new file mode 100644 index 0000000000000..8018f79479810 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_modify_singular_modes/template_tba_cfg.py @@ -0,0 +1,104 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("ppsTrackBasedAlignmentTest") + +# minimum of logs +process.MessageLogger = cms.Service("MessageLogger", + statistics = cms.untracked.vstring(), + destinations = cms.untracked.vstring('cout'), + cout = cms.untracked.PSet( + threshold = cms.untracked.string('WARNING') + ) +) + +# random seeds +process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", + ppsFastLocalSimulation = cms.PSet( + initialSeed = cms.untracked.uint32(81) + ) +) + +# data source +process.source = cms.Source("EmptySource") + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10000) +) + +# geometry +process.load("Geometry.VeryForwardGeometry.geometryRPFromDD_2018_cfi") +del(process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles[-1]) +process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles.append("CalibPPS/AlignmentRelative/test/test_with_mc/RP_Dist_Beam_Cent.xml") + +# initial alignments +process.load("CalibPPS.ESProducers.ctppsRPAlignmentCorrectionsDataESSourceXML_cfi") +process.ctppsRPAlignmentCorrectionsDataESSourceXML.MisalignedFiles = cms.vstring("CalibPPS/AlignmentRelative/test/test_modify_singular_modes/$misalignment_file") + +# simulation +process.load("CalibPPS.AlignmentRelative.ppsFastLocalSimulation_cfi") +process.ppsFastLocalSimulation.verbosity = 0 +process.ppsFastLocalSimulation.z0 = 215000 +process.ppsFastLocalSimulation.RPs = cms.vuint32(103, 104, 105, 123, 124, 125) +process.ppsFastLocalSimulation.roundToPitch = False + +# strips: pattern recognition +process.load("RecoPPS.Local.totemRPUVPatternFinder_cfi") +process.totemRPUVPatternFinder.tagRecHit = cms.InputTag("ppsFastLocalSimulation") + +# aligner +process.load("CalibPPS.AlignmentRelative.ppsStraightTrackAligner_cfi") + +process.ppsStraightTrackAligner.verbosity = 1 + +process.ppsStraightTrackAligner.tagUVPatternsStrip = cms.InputTag("totemRPUVPatternFinder") +process.ppsStraightTrackAligner.tagDiamondHits = cms.InputTag("") +process.ppsStraightTrackAligner.tagPixelHits = cms.InputTag("ppsFastLocalSimulation") + +process.ppsStraightTrackAligner.maxEvents = 1000 + +process.ppsStraightTrackAligner.rpIds = [103, 104, 105, 123, 124, 125] +process.ppsStraightTrackAligner.excludePlanes = cms.vuint32() +process.ppsStraightTrackAligner.z0 = process.ppsFastLocalSimulation.z0 + +process.ppsStraightTrackAligner.maxResidualToSigma = 100 +process.ppsStraightTrackAligner.minimumHitsPerProjectionPerRP = 3 + +process.ppsStraightTrackAligner.removeImpossible = True +process.ppsStraightTrackAligner.requireNumberOfUnits = 2 +process.ppsStraightTrackAligner.requireOverlap = False +process.ppsStraightTrackAligner.requireAtLeast3PotsInOverlap = True +process.ppsStraightTrackAligner.additionalAcceptedRPSets = "" + +process.ppsStraightTrackAligner.cutOnChiSqPerNdf = True +process.ppsStraightTrackAligner.chiSqPerNdfCut = 5000 + +process.ppsStraightTrackAligner.maxTrackAx = 1 +process.ppsStraightTrackAligner.maxTrackAy = 1 + +process.ppsStraightTrackAligner.resolveShR = True +process.ppsStraightTrackAligner.resolveShZ = False +process.ppsStraightTrackAligner.resolveRotZ = True + +process.ppsStraightTrackAligner.constraintsType = cms.string("standard") +process.ppsStraightTrackAligner.standardConstraints.units = cms.vuint32(101, 121) + +process.ppsStraightTrackAligner.algorithms = cms.vstring("Ideal", "Jan") + +process.ppsStraightTrackAligner.JanAlignmentAlgorithm.stopOnSingularModes = False + +process.ppsStraightTrackAligner.taskDataFileName = "" # results_dir + "/task_data.root" + +process.ppsStraightTrackAligner.fileNamePrefix = "$output_iteration" +process.ppsStraightTrackAligner.expandedFileNamePrefix = "" +process.ppsStraightTrackAligner.factoredFileNamePrefix = "$output_cumulative" + +process.ppsStraightTrackAligner.diagnosticsFile = "" +process.ppsStraightTrackAligner.buildDiagnosticPlots = False +process.ppsStraightTrackAligner.JanAlignmentAlgorithm.buildDiagnosticPlots = False + +# processing sequence +process.p = cms.Path( + process.ppsFastLocalSimulation + * process.totemRPUVPatternFinder + * process.ppsStraightTrackAligner +) diff --git a/CalibPPS/AlignmentRelative/test/test_with_data/RP_Dist_Beam_Cent.xml b/CalibPPS/AlignmentRelative/test/test_with_data/RP_Dist_Beam_Cent.xml new file mode 100644 index 0000000000000..2602603746637 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_data/RP_Dist_Beam_Cent.xml @@ -0,0 +1,47 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/CalibPPS/AlignmentRelative/test/test_with_data/run_test b/CalibPPS/AlignmentRelative/test/test_with_data/run_test new file mode 100755 index 0000000000000..cfc0d1d024c0a --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_data/run_test @@ -0,0 +1,4 @@ +#!/bin/bash + +./tb_submit tb_config_sector_45 -run +./tb_submit tb_config_sector_56 -run diff --git a/CalibPPS/AlignmentRelative/test/test_with_data/tb_code b/CalibPPS/AlignmentRelative/test/test_with_data/tb_code new file mode 100644 index 0000000000000..52951b8d0356f --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_data/tb_code @@ -0,0 +1,393 @@ +#!/bin/bash + +# directory to search for RECO files +reco_dir="" + +# directory with alignment code and results +top_dir=`pwd -P` + +# template for CMSSW config file +template="tb_template.py" + +# file names +cfg_file="cfg.py" +log_file="log" +job_file="job" +finished_file="finished" +time_file="timestamp" + +# settings for usin LXBATCH +bsub_options="-o /dev/null -e /dev/null" + +# name settings for subdirectory structure +subdir="tb" +label="" + +# name of geometry file +geometry="" + +# rarely changing settings for alignment module +iterationsWithOneRotZPerPot="0" +runsWithoutHorizontalRPs="" + +# alignment result that is to be used as input for the next iteration +iterationFile='results_cumulative_factored_Jan.xml' + +#---------------------------------------------------------------------------------------------------- + +job_dir="" + +function InitJobDir() +{ + job_id="$1" + job_dir="$top_dir/$job_id" + job_dir_rel="$job_id" + + rm -rf "$job_dir" + mkdir -p "$job_dir" +} + +#---------------------------------------------------------------------------------------------------- + +function CleanJobDir() +{ + echo "" + #echo "rm -rf \"$job_dir\"" +} + +#---------------------------------------------------------------------------------------------------- + +function BeginJobScript() +{ + cwd=`pwd -P` + + echo "#!/bin/sh" + echo "job_home=\"$cwd\"" + #echo "export HOME=\"$cwd\"" + echo "source \"/cvmfs/cms.cern.ch/cmsset_default.sh\"" + echo "cd \"$CMSSW_BASE\"" + echo "eval \`scram runtime -sh\`" + echo "" + echo "export CMSSW_SEARCH_PATH=\"\$CMSSW_SEARCH_PATH:$top_dir\"" + echo "export CMSSW_DATA_PATH=\"/afs\"" + echo "" + echo "job_dir=\"$job_dir\"" + echo "" + echo "cd \"\$job_home\"" +} + +#---------------------------------------------------------------------------------------------------- + +function EndJobScript() +{ + CleanJobDir + + echo "touch \"\$job_dir/$finished_file\"" +} + +#---------------------------------------------------------------------------------------------------- + +function MakeJobScript() +{ + ( + BeginJobScript + + for i in `seq $iterations` + do + echo "" + echo "### iteration $i ###" + echo "echo \"* iteration $i\"" + echo "date > \"\$job_dir/iteration$i/$time_file\"" + echo "cmsRun \"\$job_dir/iteration$i/$cfg_file\" &> \"\$job_dir/iteration$i/$log_file\"" + done + + echo "" + echo "echo \"iterations finished\"" + + EndJobScript + ) > "$job_dir/$job_file" + + chmod a+x "$job_dir/$job_file" +} + +#---------------------------------------------------------------------------------------------------- + +function AddInputFiles() +{ + run="$1" + + # query the relevant storage + if [[ "$reco_dir" == "/eos/"* ]] + then + prefix="root://eostotem.cern.ch/" + file_list_raw=`eos ls "$reco_dir" | grep "run_${run}.*.root"` + fi + + if [[ "$reco_dir" == "/castor/"* ]] + then + prefix="rfio://" + file_list_raw=`nsls "$reco_dir" | grep "run_${run}.*.root"` + fi + + if [ "$prefix" == "" ] + then + prefix="file://" + file_list_raw=`ls -1 "$reco_dir" | grep "reco_${run}.*_reco.root"` + fi + + # append the result to the file list + for file in $file_list_raw + do + if [ -n "$file_list" ] + then + file_list="$file_list,\n" + fi + + file_list="$file_list '$prefix$reco_dir/$file'" + done +} + +#---------------------------------------------------------------------------------------------------- + +function MakeCfgFiles() +{ + # build input file list + prefix="" + file_list="" + ls_list="" + + runs=($task) + for runLsSpec in "${runs[@]}" + do + local run="${runLsSpec%.*}" + local lsSpec="${runLsSpec#*.}" + + if [ -n "$ls_list" ] + then + ls_list="$ls_list, " + fi + + if [ "$lsSpec" == "$runLsSpec" ] + then + ls_list="${ls_list}\"$run:1-$run:999999\"" + else + local ls_min="${lsSpec%-*}" + local ls_max="${lsSpec#*-}" + ls_list="${ls_list}\"$run:${ls_min}-$run:${ls_max}\"" + fi + + AddInputFiles "$run" + done + + # make iteration directories + for i in `seq $iterations` + do + dir="iteration$i" + mkdir -p "$dir" + + # resolve initial alignments + alignmentFiles="" + if [ $i -eq 1 ] + then + for af in ${initial_alignment[*]} + do + #if [ ${af:0:1} != "/" ]; then af="${job_dir}/$af"; fi + + if [ -n "$alignmentFiles" ]; then alignmentFiles="$alignmentFiles, "; fi + alignmentFiles="${alignmentFiles}'${af}'" + done + else + let p=i-1 + alignmentFiles="'$job_dir_rel/iteration$p/${iterationFile}'" + fi + + # what is to be optimized in this iteration + let li=${#optimize[@]}-1 + if [ $i -le $li ]; then let li=$i-1; fi + optimize_value=${optimize[$li]} + + resolveShR=False + resolveRotZ=False + if [ -n "`echo $optimize_value|grep s`" ]; then resolveShR=True; fi + if [ -n "`echo $optimize_value|grep r`" ]; then resolveRotZ=True; fi + + # what constraints to be applied + #let li=${#useExtendedConstraints[@]}-1 + #if [ $i -le $li ]; then let li=$i-1; fi + #useExtendedConstraints_value=${useExtendedConstraints[$li]} + + # use 1 rotz per RP? + oneRotZPerPot="False" + if [ $i -le "$iterationsWithOneRotZPerPot" ]; then oneRotZPerPot="True"; fi + + # chi^2/ndf cut for this iteration + let li=${#chiSqPerNdfCut[@]}-1 + if [ $i -le $li ]; then let li=$i-1; fi + chiSqPerNdfCut_value=${chiSqPerNdfCut[$li]} + + # R/sigma cut for this iteration + let li=${#maxResidualToSigma[@]}-1 + if [ $i -le $li ]; then let li=$i-1; fi + maxResidualToSigma_value=${maxResidualToSigma[$li]} + + # max a_x cut for this iteration + let li=${#maxTrackAx[@]}-1 + if [ $i -le $li ]; then let li=$i-1; fi + maxTrackAx_value=${maxTrackAx[$li]} + + # max a_y cut for this iteration + let li=${#maxTrackAy[@]}-1 + if [ $i -le $li ]; then let li=$i-1; fi + maxTrackAy_value=${maxTrackAy[$li]} + + # build diagnotic plots? + buildDiagnosticPlots="False" + if [ "$i" -le "$iterationsWithDiagnostics" ]; then buildDiagnosticPlots="True"; fi + + # fill in template + cat "$top_dir/$template" | sed "\ + s|\$inputFiles|$file_list|;\ + s|\$lsList|$ls_list|;\ + s|\$optimize|$optimize_value|;\ + s|\$rps|$rps|;\ + s|\$excludePlanes|$excludePlanes|;\ + s|\$z0|$z0|;\ + s|\$resolveShR|$resolveShR|;\ + s|\$resolveRotZ|$resolveRotZ|;\ + s|\$constraintsType|$constraintsType|;\ + s|\$fixed_planes_shr|$fixed_planes_shr|;\ + s|\$fixed_planes_values_shr|$fixed_planes_values_shr|;\ + s|\$fixed_planes_rotz|$fixed_planes_rotz|;\ + s|\$fixed_planes_values_rotz|$fixed_planes_values_rotz|;\ + s|\$geometry|$geometry|;\ + s|\$reco_dir|$reco_dir|;\ + s|\$results_dir|$job_dir/$dir|;\ + s|\$chiSqPerNdfCut|$chiSqPerNdfCut_value|;\ + s|\$alignmentFiles|$alignmentFiles|;\ + s|\$useExtendedConstraints|$useExtendedConstraints_value|;\ + s|\$useZeroThetaRotZConstraint|$useZeroThetaRotZConstraint|;\ + s|\$minimumHitsPerProjectionPerRP|$minimumHitsPerProjectionPerRP|;\ + s|\$maxResidualToSigma|$maxResidualToSigma_value|;\ + s|\$oneRotZPerPot|$oneRotZPerPot|;\ + s|\$requireNumberOfUnits|$requireNumberOfUnits|;\ + s|\$requireOverlap|$requireOverlap|;\ + s|\$requireAtLeast3PotsInOverlap|$requireAtLeast3PotsInOverlap|;\ + s|\$buildDiagnosticPlots|$buildDiagnosticPlots|;\ + s|\$maxTrackAx|$maxTrackAx_value|;\ + s|\$maxTrackAy|$maxTrackAy_value|;\ + s|\$useEqualMeanUMeanVRotZConstraints|$useEqualMeanUMeanVRotZConstraints|;\ + s|\$final_constraints_units|$final_constraints_units|;\ + s|\$maxEvents|$maxEvents|;\ + s|\$additionalAcceptedRPSets|$additionalAcceptedRPSets|;\ + " > "$dir/$cfg_file" + done +} + +#---------------------------------------------------------------------------------------------------- + +function SubmitOrExecute() +{ + if [ "$test_only" == "y" ] + then + echo " execution inhibited (test_only=t)" + return 1 + fi + + if [ "$use_lxbatch" == "y" ] + then + #echo "bsub $job_dir/$job_file" + res=`bsub -q $queue $bsub_options "$job_dir/$job_file"` + echo " $res" + echo "$res" > "$job_dir/submitted" + else + #echo "$job_dir/$job_file" + $job_dir/$job_file & + echo " job executed in background" + fi +} + +#---------------------------------------------------------------------------------------------------- + +function AddToCaption() +{ + val="f" + if [ "$2" == "True" ]; then val="t"; fi + caption="${caption},$1=$val" +} + +#---------------------------------------------------------------------------------------------------- + +function CompileCaption() +{ + opt_lab="" + opt_prev="" + for opt in "${optimize[@]}" + do + if [ "$opt_prev" != "$opt" ] + then + if [ -n "$opt_lab" ]; then opt_lab="$opt_lab+"; fi + opt_lab="$opt_lab$opt" + fi + opt_prev="$opt" + done + + caption="ev=${maxEvents},pl=${minimumHitsPerProjectionPerRP},units=${requireNumberOfUnits}" + AddToCaption "ovlp" "$requireOverlap" + AddToCaption "3rpsInO" "$requireAtLeast3PotsInOverlap" + + caption="$caption/${opt_lab},${constraintsLabel},1rotzIt=${iterationsWithOneRotZPerPot}" + AddToCaption "eMuMvRot" "$useEqualMeanUMeanVRotZConstraints" + + caption="${caption}${label}" +} + +#---------------------------------------------------------------------------------------------------- + +function SubmitTrackBasedAlignment() +{ + CompileCaption + + for task in "${tasks[@]}" + do + runs=($task) + task_id="" + for run in "${runs[@]}" + do + if [ -n "$task_id" ] + then + task_id="$task_id," + fi + + task_id="${task_id}${run}" + done + + task_id=${task_id///} + + input_tag="$rps" + if [ -n "$excludePlanes" ] + then + input_tag="${input_tag}-excl${excludePlanes}" + fi + + id="${subdir}/${task_id}/${input_tag}/${caption}" + + echo "> $id" + + InitJobDir "$id" + cd "$job_dir" + + MakeJobScript + MakeCfgFiles + + SubmitOrExecute + done +} + +#---------------------------------------------------------------------------------------------------- + +function NewDataSet() +{ + # reset standard settings + tasks="" +} diff --git a/CalibPPS/AlignmentRelative/test/test_with_data/tb_config_sector_45 b/CalibPPS/AlignmentRelative/test/test_with_data/tb_config_sector_45 new file mode 100644 index 0000000000000..621b9f73a94c4 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_data/tb_config_sector_45 @@ -0,0 +1,61 @@ +#!/bin/bash + +# what to optimise +sector="45" +z0="-210000" +rps="3,4,5,23,24,25" + +# planes 4.4, 23.5 +excludePlanes="1981939712,2023555072" + +# event selection +minimumHitsPerProjectionPerRP=3 +requireNumberOfUnits="2" +requireOverlap="False" +requireAtLeast3PotsInOverlap="True" + +additionalAcceptedRPSets="" + +# even number limit +#maxEvents="-1" # process all events +maxEvents="1E4" + +# iteration settings +iterations=5 +optimize=("s" "s" "sr") +iterationsWithOneRotZPerPot=0 +maxResidualToSigma=(30 10 10 3) +chiSqPerNdfCut=(500 50 50 5) +maxTrackAx=("0.5E-3") +maxTrackAy=("0.5E-3") + +# constraint settings +constraintsType="standard" +constraintsLabel="std" +final_constraints_units="1,21" +useEqualMeanUMeanVRotZConstraints="True" + +# diagonistics +iterationsWithDiagnostics="1" + +# input settings +reco_dir="/eos/totem/data/offline/2018/450GeV/beta11/ZeroBias/version1" + +# geometry settings +geometry="CalibPPS/AlignmentRelative/test/test_with_data" + +initial_alignment=() + +# output settings +subdir="data" +label="" # for non-standard settings, e.g. ",test" + +NewDataSet +tasks=( + 324579.1-25 + 324579.26-50 + 324579.51-75 + 324579.76-99 +) + +SubmitTrackBasedAlignment diff --git a/CalibPPS/AlignmentRelative/test/test_with_data/tb_config_sector_56 b/CalibPPS/AlignmentRelative/test/test_with_data/tb_config_sector_56 new file mode 100644 index 0000000000000..213aa7b51e5de --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_data/tb_config_sector_56 @@ -0,0 +1,61 @@ +#!/bin/bash + +# what to optimise +sector="56" +z0="+217000" +rps="103,104,105,123,124,125" + +# nothing to exclude +excludePlanes="" + +# event selection +minimumHitsPerProjectionPerRP=3 +requireNumberOfUnits="2" +requireOverlap="False" +requireAtLeast3PotsInOverlap="True" + +additionalAcceptedRPSets="" + +# even number limit +#maxEvents="-1" # process all events +maxEvents="1E4" + +# iteration settings +iterations=5 +optimize=("s" "s" "sr") +iterationsWithOneRotZPerPot=0 +maxResidualToSigma=(30 10 10 3) +chiSqPerNdfCut=(500 50 50 5) +maxTrackAx=("0.5E-3") +maxTrackAy=("0.5E-3") + +# constraint settings +constraintsType="standard" +constraintsLabel="std" +final_constraints_units="101,121" +useEqualMeanUMeanVRotZConstraints="True" + +# diagonistics +iterationsWithDiagnostics="1" + +# input settings +reco_dir="/eos/totem/data/offline/2018/450GeV/beta11/ZeroBias/version1" + +# geometry settings +geometry="CalibPPS/AlignmentRelative/test/test_with_data" + +initial_alignment=() + +# output settings +subdir="data" +label="" # for non-standard settings, e.g. ",test" + +NewDataSet +tasks=( + 324579.1-25 + 324579.26-50 + 324579.51-75 + 324579.76-99 +) + +SubmitTrackBasedAlignment diff --git a/CalibPPS/AlignmentRelative/test/test_with_data/tb_submit b/CalibPPS/AlignmentRelative/test/test_with_data/tb_submit new file mode 100755 index 0000000000000..8486b48d440a8 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_data/tb_submit @@ -0,0 +1,97 @@ +#!/bin/bash + +#---------------------------------------------------------------------------------------------------- + +# set defaults +test_only="y" +use_lxbatch="n" +queue="1nd" + +#---------------------------------------------------------------------------------------------------- + +function PrintUsage() +{ + echo "USAGE: submit [config]" + echo "OPTIONS:" + echo " -test prepare jobs but do not execute" + echo " -run prepare jobs and execute" + echo " -local execute jobs locally" + echo " -batch, -lxbatch submit jobs to LXBATCH" + echo " -submit abbreviation of -run and -batch" + echo " -queue use LXBATCH queue , implies -batch option" +} + +#---------------------------------------------------------------------------------------------------- + +# load code +source "tb_code" + +#---------------------------------------------------------------------------------------------------- +# parse command line + +while [ -n "$1" ] +do + case "$1" in + "-h" | "--help") + PrintUsage + exit 0 + ;; + + "-test") + test_only="y" + ;; + + "-run") + test_only="n" + ;; + + "-local") + use_lxbatch="n" + ;; + + "-batch" | "-lxbatch") + use_lxbatch="y" + ;; + + "-queue") + shift + queue="$1" + use_lxbatch="y" + ;; + + "-submit") + test_only="n" + use_lxbatch="y" + ;; + + -*) + echo "ERROR: unknown option '$1'." + PrintUsage + exit 1 + ;; + + *) + if [ -n "$input_config" ] + then + echo "ERROR: only one config can be used at at time." + exit 1 + fi + input_config="$1" + ;; + + esac + + shift +done + +#---------------------------------------------------------------------------------------------------- +# source config + +if [ ! -f "$input_config" ] +then + echo "ERROR: can't read config file '$input_config'." + PrintUsage + exit 1 +fi + +source "$input_config" diff --git a/CalibPPS/AlignmentRelative/test/test_with_data/tb_template.py b/CalibPPS/AlignmentRelative/test/test_with_data/tb_template.py new file mode 100644 index 0000000000000..81f9571e4921e --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_data/tb_template.py @@ -0,0 +1,113 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.StandardSequences.Eras import eras +process = cms.Process("trackBasedAlignment", eras.Run2_2018) + +# minimum of logs +process.MessageLogger = cms.Service("MessageLogger", + statistics = cms.untracked.vstring(), + destinations = cms.untracked.vstring('cout'), + cout = cms.untracked.PSet( + threshold = cms.untracked.string('WARNING') + ) +) + +# input data +process.source = cms.Source("PoolSource", + skipBadFiles = cms.untracked.bool(True), + fileNames = cms.untracked.vstring( +$inputFiles + ), + lumisToProcess = cms.untracked.VLuminosityBlockRange($lsList), + inputCommands = cms.untracked.vstring( + "drop *", + "keep TotemRPRecHitedmDetSetVector_*_*_*", + "keep CTPPSPixelRecHitedmDetSetVector_*_*_*", + ) +) + +# geometry +process.load("Geometry.VeryForwardGeometry.geometryRPFromDD_2018_cfi") +del(process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles[-1]) +process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles.append("$geometry/RP_Dist_Beam_Cent.xml") + +# initial alignments +process.load("CalibPPS.ESProducers.ctppsRPAlignmentCorrectionsDataESSourceXML_cfi") +process.ctppsRPAlignmentCorrectionsDataESSourceXML.RealFiles = cms.vstring($alignmentFiles) + +# reco modules +process.load("RecoPPS.Local.totemRPLocalReconstruction_cff") + +process.load("RecoPPS.Local.ctppsPixelLocalReconstruction_cff") + +process.load("RecoPPS.Local.ctppsLocalTrackLiteProducer_cff") +process.ctppsLocalTrackLiteProducer.includeDiamonds = False + +# aligner +process.load("CalibPPS.AlignmentRelative.ppsStraightTrackAligner_cfi") + +process.ppsStraightTrackAligner.verbosity = 1 + +process.ppsStraightTrackAligner.tagUVPatternsStrip = cms.InputTag("totemRPUVPatternFinder") +process.ppsStraightTrackAligner.tagDiamondHits = cms.InputTag("") +process.ppsStraightTrackAligner.tagPixelHits = cms.InputTag("") +process.ppsStraightTrackAligner.tagPixelLocalTracks = cms.InputTag("ctppsPixelLocalTracks") + +process.ppsStraightTrackAligner.maxEvents = int($maxEvents) + +process.ppsStraightTrackAligner.rpIds = [$rps] +process.ppsStraightTrackAligner.excludePlanes = cms.vuint32($excludePlanes) +process.ppsStraightTrackAligner.z0 = $z0 + +process.ppsStraightTrackAligner.maxResidualToSigma = $maxResidualToSigma +process.ppsStraightTrackAligner.minimumHitsPerProjectionPerRP = $minimumHitsPerProjectionPerRP + +process.ppsStraightTrackAligner.removeImpossible = True +process.ppsStraightTrackAligner.requireNumberOfUnits = $requireNumberOfUnits +process.ppsStraightTrackAligner.requireOverlap = $requireOverlap +process.ppsStraightTrackAligner.requireAtLeast3PotsInOverlap = $requireAtLeast3PotsInOverlap +process.ppsStraightTrackAligner.additionalAcceptedRPSets = "$additionalAcceptedRPSets" + +process.ppsStraightTrackAligner.cutOnChiSqPerNdf = True +process.ppsStraightTrackAligner.chiSqPerNdfCut = $chiSqPerNdfCut + +process.ppsStraightTrackAligner.maxTrackAx = $maxTrackAx +process.ppsStraightTrackAligner.maxTrackAy = $maxTrackAy + +optimize="$optimize" +process.ppsStraightTrackAligner.resolveShR = $resolveShR +process.ppsStraightTrackAligner.resolveShZ = False +process.ppsStraightTrackAligner.resolveRotZ = $resolveRotZ + +process.ppsStraightTrackAligner.constraintsType = "standard" +process.ppsStraightTrackAligner.standardConstraints.units = cms.vuint32($final_constraints_units) +process.ppsStraightTrackAligner.oneRotZPerPot = $oneRotZPerPot +process.ppsStraightTrackAligner.useEqualMeanUMeanVRotZConstraints = $useEqualMeanUMeanVRotZConstraints + +process.ppsStraightTrackAligner.algorithms = cms.vstring("Jan") + +process.ppsStraightTrackAligner.JanAlignmentAlgorithm.stopOnSingularModes = False + +results_dir="$results_dir" + +process.ppsStraightTrackAligner.taskDataFileName = "" # results_dir + "/task_data.root" + +process.ppsStraightTrackAligner.fileNamePrefix = results_dir + "/results_iteration_" +process.ppsStraightTrackAligner.expandedFileNamePrefix = results_dir + "/results_cumulative_expanded_" +process.ppsStraightTrackAligner.factoredFileNamePrefix = results_dir + "/results_cumulative_factored_" + +process.ppsStraightTrackAligner.diagnosticsFile = results_dir + '/diagnostics.root' +process.ppsStraightTrackAligner.buildDiagnosticPlots = $buildDiagnosticPlots +process.ppsStraightTrackAligner.JanAlignmentAlgorithm.buildDiagnosticPlots = $buildDiagnosticPlots + +# processing sequence +process.p = cms.Path( + # it is important to re-run part of the reconstruction as it may influence + # the choice of rec-hits used in the alignment + process.totemRPUVPatternFinder + * process.totemRPLocalTrackFitter + * process.ctppsPixelLocalTracks + * process.ctppsLocalTrackLiteProducer + + * process.ppsStraightTrackAligner +) diff --git a/CalibPPS/AlignmentRelative/test/test_with_mc/RP_Dist_Beam_Cent.xml b/CalibPPS/AlignmentRelative/test/test_with_mc/RP_Dist_Beam_Cent.xml new file mode 100644 index 0000000000000..c3443025f2c08 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_mc/RP_Dist_Beam_Cent.xml @@ -0,0 +1,42 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/CalibPPS/AlignmentRelative/test/test_with_mc/iterations/run_test b/CalibPPS/AlignmentRelative/test/test_with_mc/iterations/run_test new file mode 100755 index 0000000000000..dbc79af54e62d --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_mc/iterations/run_test @@ -0,0 +1,28 @@ +#!/bin/bash + +iterations="4" + +mkdir -p "data" + +cd "data" + +for iter in `seq $iterations` +do + echo "* iteration $iter" + + if [ "$iter" == "1" ] + then + input_alignment="" + else + let pred=$iter-1 + input_alignment="\"CalibPPS/AlignmentRelative/test/test_with_mc/iterations/data/iteration_${pred}_cumulative_Jan.xml\"" + fi + + cat "../template_cfg.py" | sed "\ + s|\$input_alignment|$input_alignment|;\ + s|\$output_iteration|iteration_${iter}_|;\ + s|\$output_cumulative|iteration_${iter}_cumulative_|;\ + " > "iteration_${iter}_cfg.py" + + cmsRun "iteration_${iter}_cfg.py" &> "iteration_${iter}.log" +done diff --git a/CalibPPS/AlignmentRelative/test/test_with_mc/iterations/template_cfg.py b/CalibPPS/AlignmentRelative/test/test_with_mc/iterations/template_cfg.py new file mode 100644 index 0000000000000..59a9caabf4657 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_mc/iterations/template_cfg.py @@ -0,0 +1,105 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("ppsTrackBasedAlignmentTest") + +# minimum of logs +process.MessageLogger = cms.Service("MessageLogger", + statistics = cms.untracked.vstring(), + destinations = cms.untracked.vstring('cout'), + cout = cms.untracked.PSet( + threshold = cms.untracked.string('WARNING') + ) +) + +# random seeds +process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", + ppsFastLocalSimulation = cms.PSet( + initialSeed = cms.untracked.uint32(81) + ) +) + +# data source +process.source = cms.Source("EmptySource") + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10000) +) + +# geometry +process.load("Geometry.VeryForwardGeometry.geometryRPFromDD_2018_cfi") +del(process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles[-1]) +process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles.append("CalibPPS/AlignmentRelative/test/test_with_mc//RP_Dist_Beam_Cent.xml") + +# initial alignments +process.load("CalibPPS.ESProducers.ctppsRPAlignmentCorrectionsDataESSourceXML_cfi") +process.ctppsRPAlignmentCorrectionsDataESSourceXML.MisalignedFiles = cms.vstring("CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_sr.xml") +process.ctppsRPAlignmentCorrectionsDataESSourceXML.RealFiles = cms.vstring($input_alignment) + +# simulation +process.load("CalibPPS.AlignmentRelative.ppsFastLocalSimulation_cfi") +process.ppsFastLocalSimulation.verbosity = 0 +process.ppsFastLocalSimulation.z0 = 215000 +process.ppsFastLocalSimulation.RPs = cms.vuint32(103, 104, 105, 123, 124, 125) +process.ppsFastLocalSimulation.roundToPitch = False + +# strips: pattern recognition +process.load("RecoPPS.Local.totemRPUVPatternFinder_cfi") +process.totemRPUVPatternFinder.tagRecHit = cms.InputTag("ppsFastLocalSimulation") + +# aligner +process.load("CalibPPS.AlignmentRelative.ppsStraightTrackAligner_cfi") + +process.ppsStraightTrackAligner.verbosity = 1 + +process.ppsStraightTrackAligner.tagUVPatternsStrip = cms.InputTag("totemRPUVPatternFinder") +process.ppsStraightTrackAligner.tagDiamondHits = cms.InputTag("") +process.ppsStraightTrackAligner.tagPixelHits = cms.InputTag("ppsFastLocalSimulation") + +process.ppsStraightTrackAligner.maxEvents = 1000 + +process.ppsStraightTrackAligner.rpIds = [103, 104, 105, 123, 124, 125] +process.ppsStraightTrackAligner.excludePlanes = cms.vuint32() +process.ppsStraightTrackAligner.z0 = process.ppsFastLocalSimulation.z0 + +process.ppsStraightTrackAligner.maxResidualToSigma = 100 +process.ppsStraightTrackAligner.minimumHitsPerProjectionPerRP = 3 + +process.ppsStraightTrackAligner.removeImpossible = True +process.ppsStraightTrackAligner.requireNumberOfUnits = 2 +process.ppsStraightTrackAligner.requireOverlap = False +process.ppsStraightTrackAligner.requireAtLeast3PotsInOverlap = True +process.ppsStraightTrackAligner.additionalAcceptedRPSets = "" + +process.ppsStraightTrackAligner.cutOnChiSqPerNdf = True +process.ppsStraightTrackAligner.chiSqPerNdfCut = 5000 + +process.ppsStraightTrackAligner.maxTrackAx = 1 +process.ppsStraightTrackAligner.maxTrackAy = 1 + +process.ppsStraightTrackAligner.resolveShR = True +process.ppsStraightTrackAligner.resolveShZ = False +process.ppsStraightTrackAligner.resolveRotZ = True + +process.ppsStraightTrackAligner.constraintsType = cms.string("standard") +process.ppsStraightTrackAligner.standardConstraints.units = cms.vuint32(101, 121) + +process.ppsStraightTrackAligner.algorithms = cms.vstring("Ideal", "Jan") + +process.ppsStraightTrackAligner.JanAlignmentAlgorithm.stopOnSingularModes = False + +process.ppsStraightTrackAligner.taskDataFileName = "" # results_dir + "/task_data.root" + +process.ppsStraightTrackAligner.fileNamePrefix = "$output_iteration" +process.ppsStraightTrackAligner.expandedFileNamePrefix = "" +process.ppsStraightTrackAligner.factoredFileNamePrefix = "$output_cumulative" + +process.ppsStraightTrackAligner.diagnosticsFile = "" +process.ppsStraightTrackAligner.buildDiagnosticPlots = False +process.ppsStraightTrackAligner.JanAlignmentAlgorithm.buildDiagnosticPlots = False + +# processing sequence +process.p = cms.Path( + process.ppsFastLocalSimulation + * process.totemRPUVPatternFinder + * process.ppsStraightTrackAligner +) diff --git a/CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_r.xml b/CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_r.xml new file mode 100644 index 0000000000000..1f811c9e2837f --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_r.xml @@ -0,0 +1,27 @@ + + + + + + + + + + + + + + + + + + + + + + diff --git a/CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_s.xml b/CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_s.xml new file mode 100644 index 0000000000000..ddd90b2fc0c73 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_s.xml @@ -0,0 +1,27 @@ + + + + + + + + + + + + + + + + + + + + + + diff --git a/CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_sr.xml b/CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_sr.xml new file mode 100644 index 0000000000000..33eac259eafcf --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_sr.xml @@ -0,0 +1,27 @@ + + + + + + + + + + + + + + + + + + + + + + diff --git a/CalibPPS/AlignmentRelative/test/test_with_mc/simple/cfg.py b/CalibPPS/AlignmentRelative/test/test_with_mc/simple/cfg.py new file mode 100644 index 0000000000000..6806e38c73827 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_mc/simple/cfg.py @@ -0,0 +1,106 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("ppsTrackBasedAlignmentTest") + +# minimum of logs +process.MessageLogger = cms.Service("MessageLogger", + statistics = cms.untracked.vstring(), + destinations = cms.untracked.vstring('cout'), + cout = cms.untracked.PSet( + threshold = cms.untracked.string('WARNING') + ) +) + +# random seeds +process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", + ppsFastLocalSimulation = cms.PSet( + initialSeed = cms.untracked.uint32(81) + ) +) + +# data source +process.source = cms.Source("EmptySource") + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10000) +) + +# geometry +process.load("Geometry.VeryForwardGeometry.geometryRPFromDD_2018_cfi") +del(process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles[-1]) +process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles.append("CalibPPS/AlignmentRelative/test/test_with_mc//RP_Dist_Beam_Cent.xml") + +# initial alignments +process.load("CalibPPS.ESProducers.ctppsRPAlignmentCorrectionsDataESSourceXML_cfi") +process.ctppsRPAlignmentCorrectionsDataESSourceXML.MisalignedFiles = cms.vstring("CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_sr.xml") + +# simulation +process.load("CalibPPS.AlignmentRelative.ppsFastLocalSimulation_cfi") +process.ppsFastLocalSimulation.verbosity = 0 +process.ppsFastLocalSimulation.z0 = 215000 +process.ppsFastLocalSimulation.RPs = cms.vuint32(103, 104, 105, 123, 124, 125) +process.ppsFastLocalSimulation.roundToPitch = True + +# strips: pattern recognition +process.load("RecoPPS.Local.totemRPUVPatternFinder_cfi") +process.totemRPUVPatternFinder.tagRecHit = cms.InputTag("ppsFastLocalSimulation") + +# aligner +process.load("CalibPPS.AlignmentRelative.ppsStraightTrackAligner_cfi") + +process.ppsStraightTrackAligner.verbosity = 1 + +process.ppsStraightTrackAligner.tagUVPatternsStrip = cms.InputTag("totemRPUVPatternFinder") +process.ppsStraightTrackAligner.tagDiamondHits = cms.InputTag("") +process.ppsStraightTrackAligner.tagPixelHits = cms.InputTag("ppsFastLocalSimulation") + +process.ppsStraightTrackAligner.maxEvents = 1000 + +process.ppsStraightTrackAligner.rpIds = [103, 104, 105, 123, 124, 125] +process.ppsStraightTrackAligner.excludePlanes = cms.vuint32() +process.ppsStraightTrackAligner.z0 = process.ppsFastLocalSimulation.z0 + +process.ppsStraightTrackAligner.maxResidualToSigma = 100 +process.ppsStraightTrackAligner.minimumHitsPerProjectionPerRP = 3 + +process.ppsStraightTrackAligner.removeImpossible = True +process.ppsStraightTrackAligner.requireNumberOfUnits = 2 +process.ppsStraightTrackAligner.requireOverlap = False +process.ppsStraightTrackAligner.requireAtLeast3PotsInOverlap = True +process.ppsStraightTrackAligner.additionalAcceptedRPSets = "" + +process.ppsStraightTrackAligner.cutOnChiSqPerNdf = True +process.ppsStraightTrackAligner.chiSqPerNdfCut = 5000 + +process.ppsStraightTrackAligner.maxTrackAx = 1 +process.ppsStraightTrackAligner.maxTrackAy = 1 + +process.ppsStraightTrackAligner.resolveShR = True +process.ppsStraightTrackAligner.resolveShZ = False +process.ppsStraightTrackAligner.resolveRotZ = True + +process.ppsStraightTrackAligner.constraintsType = cms.string("standard") +process.ppsStraightTrackAligner.standardConstraints.units = cms.vuint32(101, 121) +process.ppsStraightTrackAligner.useEqualMeanUMeanVRotZConstraints = False +process.ppsStraightTrackAligner.oneRotZPerPot = True + +process.ppsStraightTrackAligner.algorithms = cms.vstring("Ideal", "Jan") + +process.ppsStraightTrackAligner.JanAlignmentAlgorithm.stopOnSingularModes = False + +process.ppsStraightTrackAligner.taskDataFileName = "" # results_dir + "/task_data.root" + +process.ppsStraightTrackAligner.fileNamePrefix = "results_" +process.ppsStraightTrackAligner.expandedFileNamePrefix = "results_cumulative_expanded_" +process.ppsStraightTrackAligner.factoredFileNamePrefix = "results_cumulative_factored_" + +process.ppsStraightTrackAligner.diagnosticsFile = 'diagnostics.root' +process.ppsStraightTrackAligner.buildDiagnosticPlots = True +process.ppsStraightTrackAligner.JanAlignmentAlgorithm.buildDiagnosticPlots = True + +# processing sequence +process.p = cms.Path( + process.ppsFastLocalSimulation + * process.totemRPUVPatternFinder + * process.ppsStraightTrackAligner +) diff --git a/CalibPPS/AlignmentRelative/test/test_with_mc/simple/run_test b/CalibPPS/AlignmentRelative/test/test_with_mc/simple/run_test new file mode 100755 index 0000000000000..398b31e8f4264 --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_mc/simple/run_test @@ -0,0 +1,7 @@ +#!/bin/bash + +mkdir -p "data" + +cd "data" + +cmsRun "../cfg.py" &> "log" diff --git a/CalibPPS/AlignmentRelative/test/test_with_mc/statistics/run_test b/CalibPPS/AlignmentRelative/test/test_with_mc/statistics/run_test new file mode 100755 index 0000000000000..7555643fc57eb --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_mc/statistics/run_test @@ -0,0 +1,47 @@ +#!/bin/bash + +events=( + "1E2" + "4E2" + "1E3" + "4E3" + "1E5" +) + +repetitions="10" + +#---------------------------------------------------------------------------------------------------- + +for n_events in ${events[*]} +do + for rep in `seq $repetitions` + do + dir="data/sector:56/opt:sr/trdist:default/$n_events/$rep/iteration:1" + + echo "* $dir" + + mkdir -p "$dir" + + cat "template_cfg.py" | sed "\ + s|\$seed|$rep|;\ + s|\$max_events|$n_events|;\ + s|\$output_iteration|results_iteration_|;\ + s|\$output_expanded|results_expanded_|;\ + s|\$output_factored|results_factored_|;\ + " > "$dir/cfg.py" + + cd "$dir" + cmsRun "cfg.py" &> "log" & + cd - &> /dev/null + done +done + +wait + +echo "* analysis" +cd "data" +$CMSSW_BASE/test/$SCRAM_ARCH/ppsTrackBasedAligAnalRes \ + ./results_factored_Jan.xml\ + ./results_factored_Ideal.xml\ + ./results_expanded_Jan.xml\ + ./results_expanded_Ideal.xml diff --git a/CalibPPS/AlignmentRelative/test/test_with_mc/statistics/template_cfg.py b/CalibPPS/AlignmentRelative/test/test_with_mc/statistics/template_cfg.py new file mode 100644 index 0000000000000..0defa0695479e --- /dev/null +++ b/CalibPPS/AlignmentRelative/test/test_with_mc/statistics/template_cfg.py @@ -0,0 +1,106 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("ppsTrackBasedAlignmentTest") + +# minimum of logs +process.MessageLogger = cms.Service("MessageLogger", + statistics = cms.untracked.vstring(), + destinations = cms.untracked.vstring('cout'), + cout = cms.untracked.PSet( + threshold = cms.untracked.string('WARNING') + ) +) + +# random seeds +process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", + ppsFastLocalSimulation = cms.PSet( + initialSeed = cms.untracked.uint32($seed) + ) +) + +# data source +process.source = cms.Source("EmptySource") + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10000000) +) + +# geometry +process.load("Geometry.VeryForwardGeometry.geometryRPFromDD_2018_cfi") +del(process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles[-1]) +process.XMLIdealGeometryESSource_CTPPS.geomXMLFiles.append("CalibPPS/AlignmentRelative/test/test_with_mc//RP_Dist_Beam_Cent.xml") + +# initial alignments +process.load("CalibPPS.ESProducers.ctppsRPAlignmentCorrectionsDataESSourceXML_cfi") +process.ctppsRPAlignmentCorrectionsDataESSourceXML.MisalignedFiles = cms.vstring("CalibPPS/AlignmentRelative/test/test_with_mc/mis_alignment_sr.xml") + +# simulation +process.load("CalibPPS.AlignmentRelative.ppsFastLocalSimulation_cfi") +process.ppsFastLocalSimulation.verbosity = 0 +process.ppsFastLocalSimulation.z0 = 215000 +process.ppsFastLocalSimulation.RPs = cms.vuint32(103, 104, 105, 123, 124, 125) +process.ppsFastLocalSimulation.roundToPitch = True + +# strips: pattern recognition +process.load("RecoPPS.Local.totemRPUVPatternFinder_cfi") +process.totemRPUVPatternFinder.tagRecHit = cms.InputTag("ppsFastLocalSimulation") + +# aligner +process.load("CalibPPS.AlignmentRelative.ppsStraightTrackAligner_cfi") + +process.ppsStraightTrackAligner.verbosity = 1 + +process.ppsStraightTrackAligner.tagUVPatternsStrip = cms.InputTag("totemRPUVPatternFinder") +process.ppsStraightTrackAligner.tagDiamondHits = cms.InputTag("") +process.ppsStraightTrackAligner.tagPixelHits = cms.InputTag("ppsFastLocalSimulation") + +process.ppsStraightTrackAligner.maxEvents = int($max_events) + +process.ppsStraightTrackAligner.rpIds = [103, 104, 105, 123, 124, 125] +process.ppsStraightTrackAligner.excludePlanes = cms.vuint32() +process.ppsStraightTrackAligner.z0 = process.ppsFastLocalSimulation.z0 + +process.ppsStraightTrackAligner.maxResidualToSigma = 100 +process.ppsStraightTrackAligner.minimumHitsPerProjectionPerRP = 3 + +process.ppsStraightTrackAligner.removeImpossible = True +process.ppsStraightTrackAligner.requireNumberOfUnits = 2 +process.ppsStraightTrackAligner.requireOverlap = False +process.ppsStraightTrackAligner.requireAtLeast3PotsInOverlap = True +process.ppsStraightTrackAligner.additionalAcceptedRPSets = "" + +process.ppsStraightTrackAligner.cutOnChiSqPerNdf = True +process.ppsStraightTrackAligner.chiSqPerNdfCut = 5000 + +process.ppsStraightTrackAligner.maxTrackAx = 1 +process.ppsStraightTrackAligner.maxTrackAy = 1 + +process.ppsStraightTrackAligner.resolveShR = True +process.ppsStraightTrackAligner.resolveShZ = False +process.ppsStraightTrackAligner.resolveRotZ = True + +process.ppsStraightTrackAligner.constraintsType = cms.string("standard") +process.ppsStraightTrackAligner.standardConstraints.units = cms.vuint32(101, 121) + +process.ppsStraightTrackAligner.algorithms = cms.vstring("Ideal", "Jan") + +process.ppsStraightTrackAligner.JanAlignmentAlgorithm.stopOnSingularModes = False + +process.ppsStraightTrackAligner.taskDataFileName = "" + +process.ppsStraightTrackAligner.fileNamePrefix = "$output_iteration" +process.ppsStraightTrackAligner.expandedFileNamePrefix = "$output_expanded" +process.ppsStraightTrackAligner.factoredFileNamePrefix = "$output_factored" + +process.ppsStraightTrackAligner.saveXMLUncertainties = True + +process.ppsStraightTrackAligner.diagnosticsFile = "" +process.ppsStraightTrackAligner.buildDiagnosticPlots = False +process.ppsStraightTrackAligner.JanAlignmentAlgorithm.buildDiagnosticPlots = False + +# processing sequence +process.p = cms.Path( + process.ppsFastLocalSimulation + * process.totemRPUVPatternFinder + * process.ppsStraightTrackAligner +) diff --git a/CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionData.h b/CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionData.h index 42000b087e343..6f9283888f39e 100644 --- a/CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionData.h +++ b/CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionData.h @@ -69,6 +69,20 @@ class CTPPSRPAlignmentCorrectionData { double rot_x_unc, rot_y_unc, rot_z_unc; public: + CTPPSRPAlignmentCorrectionData() + : sh_x(0.), + sh_y(0.), + sh_z(0.), + sh_x_unc(0.), + sh_y_unc(0.), + sh_z_unc(0.), + rot_x(0.), + rot_y(0.), + rot_z(0.), + rot_x_unc(0.), + rot_y_unc(0.), + rot_z_unc(0.) {} + /// full constructor, shifts in mm, rotations in rad CTPPSRPAlignmentCorrectionData(double _sh_x, double _sh_x_u, @@ -84,12 +98,7 @@ class CTPPSRPAlignmentCorrectionData { double _rot_z_u); /// no uncertainty constructor, shifts in mm, rotation in rad - CTPPSRPAlignmentCorrectionData(double _sh_x = 0., - double _sh_y = 0., - double _sh_z = 0., - double _rot_x = 0., - double _rot_y = 0., - double rot_z = 0.); + CTPPSRPAlignmentCorrectionData(double _sh_x, double _sh_y, double _sh_z, double _rot_x, double _rot_y, double rot_z); inline double getShX() const { return sh_x; } inline void setShX(const double &v) { sh_x = v; } @@ -136,9 +145,7 @@ class CTPPSRPAlignmentCorrectionData { } /// merges (cumulates) alignements - /// match between x, y and read-out shifts is not checked - /// \param sumErrors if it is true, old and new alignment uncertainties are summed (in quadrature) - /// if it is false, the uncertainties of the parameter (i.e. not the object) will be used + /// \param sumErrors if true, uncertainties are summed in quadrature, otherwise the uncertainties of this are not changed /// With the add... switches one can control which corrections are added. void add(const CTPPSRPAlignmentCorrectionData &, bool sumErrors = true, bool addSh = true, bool addRot = true); diff --git a/CondFormats/PPSObjects/src/CTPPSRPAlignmentCorrectionsData.cc b/CondFormats/PPSObjects/src/CTPPSRPAlignmentCorrectionsData.cc index df3c203ba6e02..79659a9032d16 100644 --- a/CondFormats/PPSObjects/src/CTPPSRPAlignmentCorrectionsData.cc +++ b/CondFormats/PPSObjects/src/CTPPSRPAlignmentCorrectionsData.cc @@ -47,19 +47,18 @@ CTPPSRPAlignmentCorrectionData CTPPSRPAlignmentCorrectionsData::getSensorCorrect CTPPSRPAlignmentCorrectionData CTPPSRPAlignmentCorrectionsData::getFullSensorCorrection(unsigned int id, bool useRPErrors) const { + // by default empty correction CTPPSRPAlignmentCorrectionData align_corr; - // try to get alignment correction of the full RP - auto rpIt = rps_.find(CTPPSDetId(id).rpId()); - if (rpIt != rps_.end()) - align_corr = rpIt->second; - - // try to get sensor alignment correction + // if found, add sensor correction (with its uncertainty) auto sIt = sensors_.find(id); - - // merge the corrections if (sIt != sensors_.end()) - align_corr.add(sIt->second, useRPErrors); + align_corr.add(sIt->second, true); + + // if found, add RP correction (depending on the flag, with or without its uncertainty) + auto rpIt = rps_.find(CTPPSDetId(id).rpId()); + if (rpIt != rps_.end()) + align_corr.add(rpIt->second, useRPErrors); return align_corr; }