MillePedeAlignmentAlgorithm/test/AlignmentRcdChecker.cpp

0001 // \file AlignmentRcdChecker.cpp
0002 //
0003 // \author    : Marco Musich
0004 // Revision   : $Revision: 1.1 $
0005 // last update: $Date: 2022/05/11 14:44:00 $
0006 // by         : $Author: musich $
0007
0008 // system includes
0009 #include <string>
0010 #include <map>
0011 #include <ranges>
0012 #include <vector>
0013 #include <tuple>
0014 #include <iterator>
0015
0016 // user includes
0017 #include "CondFormats/Alignment/interface/AlignTransform.h"
0018 #include "CondFormats/Alignment/interface/Alignments.h"
0019 #include "CondFormats/AlignmentRecord/interface/CSCAlignmentRcd.h"
0020 #include "CondFormats/AlignmentRecord/interface/DTAlignmentRcd.h"
0021 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
0022 #include "DataFormats/DetId/interface/DetId.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024 #include "FWCore/Framework/interface/ESWatcher.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0029
0030 // ROOT includes
0031 #include <TMath.h>
0032
0033 namespace checker {
0034   template <typename T1, typename T2>
0035   auto zip(const std::vector<T1>& v1, const std::vector<T2>& v2) {
0036     std::vector<std::tuple<T1, T2>> result;
0037
0038     auto minSize = std::min(v1.size(), v2.size());
0039     for (size_t i = 0; i < minSize; ++i) {
0040       result.emplace_back(v1[i], v2[i]);
0041     }
0042
0043     return result;
0044   }
0045 }  // namespace checker
0046
0047 class AlignmentRcdChecker : public edm::one::EDAnalyzer<> {
0048 public:
0049   explicit AlignmentRcdChecker(const edm::ParameterSet& iConfig);
0050   ~AlignmentRcdChecker() = default;
0051
0052   virtual void analyze(const edm::Event& evt, const edm::EventSetup& evtSetup);
0053
0054 private:
0055   void inspectRecord(const std::string& rcdname,
0056                      const edm::Event& evt,
0057                      const Alignments* refAlignments,
0058                      const Alignments* alignments);
0059
0060   const bool verbose_;
0061   const bool compareStrict_;
0062   const std::string label_;
0063
0064   const edm::ESGetToken<Alignments, TrackerAlignmentRcd> tkAliTokenRef_;
0065   const edm::ESGetToken<Alignments, TrackerAlignmentRcd> tkAliTokenNew_;
0066
0067   static constexpr double strictTolerance_ = 1e-4;  // if in cm (i.e. 1 micron)
0068 };
0069
0070 AlignmentRcdChecker::AlignmentRcdChecker(const edm::ParameterSet& iConfig)
0071     : verbose_(iConfig.getParameter<bool>("verbose")),
0072       compareStrict_(iConfig.getParameter<bool>("compareStrict")),
0073       label_(iConfig.getParameter<std::string>("label")),
0074       tkAliTokenRef_(esConsumes()),
0075       tkAliTokenNew_(esConsumes(edm::ESInputTag{"", label_})) {}
0076
0077 void AlignmentRcdChecker::analyze(const edm::Event& evt, const edm::EventSetup& evtSetup) {
0078   const Alignments* alignmentsRef = &evtSetup.getData(tkAliTokenRef_);
0079   const Alignments* alignmentsNew = &evtSetup.getData(tkAliTokenNew_);
0080   inspectRecord("TrackerAlignmentRcd", evt, alignmentsRef, alignmentsNew);
0081 }
0082
0083 void AlignmentRcdChecker::inspectRecord(const std::string& rcdname,
0084                                         const edm::Event& evt,
0085                                         const Alignments* refAlignments,
0086                                         const Alignments* alignments) {
0087   edm::LogPrint("inspectRecord") << rcdname << " content starting from run " << evt.run();
0088   edm::LogPrint("inspectRecord") << " with " << alignments->m_align.size() << " entries";
0089
0090   if (refAlignments && alignments) {
0091     double meanX = 0, rmsX = 0;
0092     double meanY = 0, rmsY = 0;
0093     double meanZ = 0, rmsZ = 0;
0094     double meanR = 0, rmsR = 0;
0095     double dPhi, meanPhi = 0, rmsPhi = 0;
0096
0097     for (const auto& [alignment, refAlignment] : checker::zip(alignments->m_align, refAlignments->m_align)) {
0098       auto delta = alignment.translation() - refAlignment.translation();
0099
0100       meanX += delta.x();
0101       rmsX += pow(delta.x(), 2);
0102
0103       meanY += delta.y();
0104       rmsY += pow(delta.y(), 2);
0105
0106       meanZ += delta.z();
0107       rmsZ += pow(delta.z(), 2);
0108
0109       meanR += delta.perp();
0110       rmsR += pow(delta.perp(), 2);
0111
0112       dPhi = alignment.translation().phi() - refAlignment.translation().phi();
0113       if (dPhi > M_PI)
0114         dPhi -= 2.0 * M_PI;
0115       if (dPhi < -M_PI)
0116         dPhi += 2.0 * M_PI;
0117
0118       meanPhi += dPhi;
0119       rmsPhi += std::pow(dPhi, 2);
0120     }
0121
0122     meanX /= alignments->m_align.size();
0123     rmsX /= alignments->m_align.size();
0124     meanY /= alignments->m_align.size();
0125     rmsY /= alignments->m_align.size();
0126     meanZ /= alignments->m_align.size();
0127     rmsZ /= alignments->m_align.size();
0128     meanR /= alignments->m_align.size();
0129     rmsR /= alignments->m_align.size();
0130     meanPhi /= alignments->m_align.size();
0131     rmsPhi /= alignments->m_align.size();
0132
0133     if (verbose_) {
0134       edm::LogPrint("inspectRecord") << "  Compared to previous record:";
0135       edm::LogPrint("inspectRecord") << "    mean X shift:   " << std::setw(12) << std::scientific
0136                                      << std::setprecision(3) << meanX << "  [cm] (RMS = " << sqrt(rmsX) << " [cm])";
0137       edm::LogPrint("inspectRecord") << "    mean Y shift:   " << std::setw(12) << std::scientific
0138                                      << std::setprecision(3) << meanY << "  [cm] (RMS = " << sqrt(rmsY) << " [cm])";
0139       edm::LogPrint("inspectRecord") << "    mean Z shift:   " << std::setw(12) << std::scientific
0140                                      << std::setprecision(3) << meanZ << "  [cm] (RMS = " << sqrt(rmsZ) << " [cm])";
0141       edm::LogPrint("inspectRecord") << "    mean R shift:   " << std::setw(12) << std::scientific
0142                                      << std::setprecision(3) << meanR << "  [cm] (RMS = " << sqrt(rmsR) << " [cm])";
0143       edm::LogPrint("inspectRecord") << "    mean Phi shift: " << std::setw(12) << std::scientific
0144                                      << std::setprecision(3) << meanPhi << " [rad] (RMS = " << sqrt(rmsPhi)
0145                                      << " [rad])";
0146     }  // verbose
0147
0148     if (compareStrict_) {
0149       // do not let any of the coordinates to fluctuate less then 1um
0150       assert(meanX < strictTolerance_);    // 1 micron
0151       assert(meanY < strictTolerance_);    // 1 micron
0152       assert(meanZ < strictTolerance_);    // 1 micron
0153       assert(meanR < strictTolerance_);    // 1 micron
0154       assert(meanPhi < strictTolerance_);  // 10 micro-rad
0155     }
0156
0157   } else {
0158     throw cms::Exception("Missing Input Data") << "Could not retrieve the input alignments to compare \n";
0159   }
0160 }
0161
0162 //define this as a plug-in
0163 DEFINE_FWK_MODULE(AlignmentRcdChecker);