BeamSpotProducer/plugins/BeamSpotCompatibilityChecker.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoVertex/BeamSpotProducer
0004 // Class:      BeamSpotCompatibilityChecker
0005 //
0006 /**\class BeamSpotCompatibilityChecker BeamSpotCompatibilityChecker.cc RecoVertex/BeamSpotProducer/plugins/BeamSpotCompatibilityChecker.cc
0007
0008  Description: Class to check the compatibility between the BeamSpot payload in the database and the one in the file
0009
0010  Implementation:
0011      Makes use of the Significance struct to establish how compatible are the data members of the two BeamSpots in input
0012 */
0013 //
0014 // Original Author:  Marco Musich
0015 //         Created:  Thu, 23 Apr 2020 09:00:45 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "CondFormats/BeamSpotObjects/interface/BeamSpotObjects.h"
0024 #include "CondFormats/DataRecord/interface/BeamSpotObjectsRcd.h"
0025 #include "CondFormats/DataRecord/interface/BeamSpotTransientObjectsRcd.h"
0026 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0027 #include "DataFormats/Common/interface/Handle.h"
0028 #include "DataFormats/Scalers/interface/BeamSpotOnline.h"
0029 #include "FWCore/Framework/interface/ESHandle.h"
0030 #include "FWCore/Framework/interface/ESWatcher.h"
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "FWCore/Framework/interface/EventSetup.h"
0033 #include "FWCore/Framework/interface/Frameworkfwd.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 #include "FWCore/Framework/interface/global/EDAnalyzer.h"
0036 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0037 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0040 #include "FWCore/Utilities/interface/ESGetToken.h"
0041 #include "FWCore/Utilities/interface/InputTag.h"
0042 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0043 //
0044 // class declaration
0045 //
0046
0047 // ancillary struct to check compatibility
0048 namespace {
0049
0050   struct Significance {
0051     Significance(const double& a, const double& b, const double& errA, const double& errB, const std::string& var)
0052         : m_A(a), m_B(b), m_ErrA(errA), m_ErrB(errB), m_var(var) {
0053       if (m_ErrA == 0 && m_ErrB == 0) {
0054         edm::LogError("LogicalError") << "Can't calculate significance when both errors are zero!" << std::endl;
0055       }
0056       m_combinedError = std::sqrt((m_ErrA * m_ErrA) + (m_ErrB * m_ErrB));
0057     }
0058
0059     inline void printBeamSpotComparison(
0060         const std::string& varName, double A, double ErrA, double B, double ErrB, double combinedError) {
0061       edm::LogPrint("BeamSpotCompatibilityChecker")
0062           << std::fixed << std::setprecision(6) << std::left << std::setw(10) << " " + varName << std::right
0063           << std::setw(12) << A << std::setw(14) << ErrA << std::setw(12) << B << std::setw(14) << ErrB << std::setw(14)
0064           << std::abs(A - B) << std::setw(14) << combinedError << std::setw(12) << std::abs(A - B) / combinedError
0065           << std::endl;
0066     }
0067
0068     float getSig(const bool verbose) {
0069       if (verbose) {
0070         printBeamSpotComparison(m_var, m_A, m_ErrA, m_B, m_ErrB, m_combinedError);
0071       }
0072       return std::abs(m_A - m_B) / m_combinedError;
0073     }
0074
0075   private:
0076     double m_A;
0077     double m_B;
0078     double m_ErrA;
0079     double m_ErrB;
0080     double m_ErrAB;
0081     std::string m_var;
0082     double m_combinedError;
0083   };
0084 }  // namespace
0085
0086 class BeamSpotCompatibilityChecker : public edm::global::EDAnalyzer<> {
0087 public:
0088   explicit BeamSpotCompatibilityChecker(const edm::ParameterSet&);
0089   ~BeamSpotCompatibilityChecker() override = default;
0090
0091   void analyze(edm::StreamID, edm::Event const&, edm::EventSetup const&) const final;
0092   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0093   static std::array<float, 6> compareBS(const reco::BeamSpot& BSA, const reco::BeamSpot& BSB, const bool verbose);
0094   static double computeBeamSpotCompatibility(const reco::BeamSpot& beamSpot1, const reco::BeamSpot& beamSpot2);
0095
0096 private:
0097   // ----------member data ---------------------------
0098   static constexpr int cmToum = 10000;
0099   const double warningThreshold_;
0100   const double throwingThreshold_;
0101   const bool verbose_;
0102   const bool dbFromEvent_;
0103
0104   // beamSpot from the file
0105   const edm::EDGetTokenT<reco::BeamSpot> bsFromFileToken_;
0106
0107   // switch to decide with record to take
0108   bool useTransientRecord_;
0109
0110   // beamSpot from the DB (object record)
0111   edm::ESGetToken<BeamSpotObjects, BeamSpotObjectsRcd> bsFromDBToken_;
0112
0113   // beamspot from the DB (transient record)
0114   edm::ESGetToken<BeamSpotObjects, BeamSpotTransientObjectsRcd> bsFromTransientDBToken_;
0115
0116   // beamSpot from the DB (via event)
0117   edm::EDGetTokenT<reco::BeamSpot> dbBSfromEventToken_;
0118
0119   template <typename RecordT>
0120   reco::BeamSpot getBeamSpotFromES(edm::EventSetup const& iSetup,
0121                                    edm::ESGetToken<BeamSpotObjects, RecordT> const& token) const {
0122     edm::ESHandle<BeamSpotObjects> beamhandle = iSetup.getHandle(token);
0123     const BeamSpotObjects* aSpot = beamhandle.product();
0124
0125     // translate from BeamSpotObjects to reco::BeamSpot
0126     reco::BeamSpot::Point apoint(aSpot->x(), aSpot->y(), aSpot->z());
0127
0128     reco::BeamSpot::CovarianceMatrix matrix;
0129     for (int i = 0; i < reco::BeamSpot::dimension; ++i)
0130       for (int j = 0; j < reco::BeamSpot::dimension; ++j)
0131         matrix(i, j) = aSpot->covariance(i, j);
0132
0133     // this assume beam width same in x and y
0134     return reco::BeamSpot(apoint, aSpot->sigmaZ(), aSpot->dxdz(), aSpot->dydz(), aSpot->beamWidthX(), matrix);
0135   }
0136 };
0137
0138 //
0139 // constructors and destructor
0140 //
0141 BeamSpotCompatibilityChecker::BeamSpotCompatibilityChecker(const edm::ParameterSet& iConfig)
0142     : warningThreshold_(iConfig.getParameter<double>("warningThr")),
0143       throwingThreshold_(iConfig.getParameter<double>("errorThr")),
0144       verbose_(iConfig.getUntrackedParameter<bool>("verbose", false)),
0145       dbFromEvent_(iConfig.getParameter<bool>("dbFromEvent")),
0146       bsFromFileToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("bsFromFile"))) {
0147   //now do what ever initialization is needed
0148   if (warningThreshold_ > throwingThreshold_) {
0149     throw cms::Exception("ConfigurationError")
0150         << __PRETTY_FUNCTION__ << "\n Warning threshold (" << warningThreshold_
0151         << ") cannot be smaller than the throwing threshold (" << throwingThreshold_ << ")" << std::endl;
0152   }
0153
0154   if (dbFromEvent_) {
0155     edm::LogInfo("BeamSpotCompatibilityChecker")
0156         << "The Database Beam Spot is going to be taken from the File via BeamSpotProducer!";
0157     dbBSfromEventToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("bsFromDB"));
0158   } else {
0159     useTransientRecord_ = iConfig.getParameter<bool>("useTransientRecord");
0160     if (useTransientRecord_) {
0161       edm::LogInfo("BeamSpotCompatibilityChecker") << "Using BeamSpot from BeamSpotTransientObjectsRcd.";
0162       bsFromTransientDBToken_ = esConsumes<BeamSpotObjects, BeamSpotTransientObjectsRcd>();
0163     } else {
0164       edm::LogInfo("BeamSpotCompatibilityChecker") << "Using BeamSpot from BeamSpotObjectsRcd.";
0165       bsFromDBToken_ = esConsumes<BeamSpotObjects, BeamSpotObjectsRcd>();
0166     }
0167     dbBSfromEventToken_ = edm::EDGetTokenT<reco::BeamSpot>();
0168   }
0169 }
0170
0171 //
0172 // member functions
0173 //
0174
0175 // ------------ method called for each event  ------------
0176 void BeamSpotCompatibilityChecker::analyze(edm::StreamID sid,
0177                                            edm::Event const& iEvent,
0178                                            edm::EventSetup const& iSetup) const {
0179   using namespace edm;
0180   reco::BeamSpot spotFile, spotDB;
0181
0182   edm::Handle<reco::BeamSpot> beamSpotFromFileHandle;
0183   iEvent.getByToken(bsFromFileToken_, beamSpotFromFileHandle);
0184   spotFile = *beamSpotFromFileHandle;
0185
0186   double file_BSx0 = spotFile.x0();
0187   double file_BSy0 = spotFile.y0();
0188   double file_BSz0 = spotFile.z0();
0189   double file_Beamsigmaz = spotFile.sigmaZ();
0190   double file_BeamWidthX = spotFile.BeamWidthX();
0191   double file_BeamWidthY = spotFile.BeamWidthY();
0192
0193   if (!dbFromEvent_) {
0194     if (useTransientRecord_) {
0195       spotDB = getBeamSpotFromES(iSetup, bsFromTransientDBToken_);
0196     } else {
0197       spotDB = getBeamSpotFromES(iSetup, bsFromDBToken_);
0198     }
0199   } else {
0200     // take the DB beamspot from the event (different label)
0201     edm::Handle<reco::BeamSpot> beamSpotFromDBHandle;
0202     iEvent.getByToken(dbBSfromEventToken_, beamSpotFromDBHandle);
0203     spotDB = *beamSpotFromDBHandle;
0204   }
0205
0206   double db_BSx0 = spotDB.x0();
0207   double db_BSy0 = spotDB.y0();
0208   double db_BSz0 = spotDB.z0();
0209   double db_Beamsigmaz = spotDB.sigmaZ();
0210   double db_BeamWidthX = spotDB.BeamWidthX();
0211   double db_BeamWidthY = spotDB.BeamWidthY();
0212
0213   double deltaX0 = file_BSx0 - db_BSx0;
0214   double deltaY0 = file_BSy0 - db_BSy0;
0215   double deltaZ0 = file_BSz0 - db_BSz0;
0216   double deltaSigmaX = file_BeamWidthX - db_BeamWidthX;
0217   double deltaSigmaY = file_BeamWidthY - db_BeamWidthY;
0218   double deltaSigmaZ = file_Beamsigmaz - db_Beamsigmaz;
0219
0220   if (verbose_) {
0221     edm::LogPrint("BeamSpotCompatibilityChecker") << "BS from DB:   \n" << spotDB << std::endl;
0222     edm::LogPrint("BeamSpotCompatibilityChecker") << "BS from File: \n" << spotFile << std::endl;
0223   }
0224
0225   auto significances = compareBS(spotDB, spotFile, verbose_);
0226   std::vector<std::string> labels = {"x0", "y0", "z0", "sigmaX", "sigmaY", "sigmaZ"};
0227
0228   std::string msg = " |delta X0|     = " + std::to_string(std::abs(deltaX0) * cmToum) +
0229                     " um\n |delta Y0|     = " + std::to_string(std::abs(deltaY0) * cmToum) +
0230                     " um\n |delta Z0|     = " + std::to_string(std::abs(deltaZ0)) +
0231                     " cm\n |delta sigmaX| = " + std::to_string(std::abs(deltaSigmaX) * cmToum) +
0232                     " um\n |delta sigmaY| = " + std::to_string(std::abs(deltaSigmaY) * cmToum) +
0233                     " um\n |delta sigmaZ| = " + std::to_string(std::abs(deltaSigmaZ)) + " cm";
0234   /*
0235   if (verbose_) {
0236     edm::LogPrint("BeamSpotCompatibilityChecker") << msg.c_str() << std::endl;
0237   }
0238   */
0239
0240   for (unsigned int i = 0; i < 3; i++) {
0241     auto sig = significances.at(i);
0242     if (sig > throwingThreshold_) {
0243       edm::LogError("BeamSpotCompatibilityChecker") << msg.c_str() << std::endl;
0244       throw cms::Exception("BeamSpotCompatibilityChecker")
0245           << "   [" << __PRETTY_FUNCTION__ << "]\n   DB-File BeamSpot " << labels.at(i) << " distance significance is "
0246           << sig << ", exceeds the threshold of " << throwingThreshold_ << "!" << std::endl;
0247     } else if (sig > warningThreshold_) {
0248       edm::LogWarning("BeamSpotCompatibilityChecker") << msg.c_str() << std::endl;
0249       edm::LogWarning("BeamSpotCompatibilityChecker")
0250           << "   [" << __PRETTY_FUNCTION__ << "]\n   DB-File BeamSpot " << labels.at(i) << " distance significance is "
0251           << sig << ", exceeds the threshold of " << warningThreshold_ << "!" << std::endl;
0252     }
0253   }
0254 }
0255
0256 std::array<float, 6> BeamSpotCompatibilityChecker::compareBS(const reco::BeamSpot& spotA,
0257                                                              const reco::BeamSpot& spotB,
0258                                                              const bool verbose) {
0259   if (verbose) {
0260     edm::LogPrint("BeamSpotCompatibilityChecker")
0261         << std::fixed << std::setprecision(6) << std::left << std::setw(10) << " Var" << std::right << std::setw(12)
0262         << "DB" << std::setw(14) << "+/-" << std::setw(12) << "File" << std::setw(14) << "+/-" << std::setw(14)
0263         << "|delta|" << std::setw(14) << "+/-" << std::setw(12) << "Sig";
0264     edm::LogPrint("BeamSpotCompatibilityChecker") << std::string(102, '-');
0265   }
0266
0267   // Lambda to calculate the significance
0268   auto calcSignificance = [&](auto a, auto b, auto aErr, auto bErr, auto var) {
0269     return Significance(a, b, aErr, bErr, var).getSig(verbose);
0270   };
0271
0272   // Populate the array using the lambda
0273   std::array<float, 6> ret = {
0274       {calcSignificance(spotA.x0(), spotB.x0(), spotA.x0Error(), spotB.x0Error(), "x"),
0275        calcSignificance(spotA.y0(), spotB.y0(), spotA.y0Error(), spotB.y0Error(), "y"),
0276        calcSignificance(spotA.z0(), spotB.z0(), spotA.z0Error(), spotB.z0Error(), "z"),
0277        calcSignificance(
0278            spotA.BeamWidthX(), spotB.BeamWidthX(), spotA.BeamWidthXError(), spotB.BeamWidthXError(), "widthX"),
0279        calcSignificance(
0280            spotA.BeamWidthY(), spotB.BeamWidthY(), spotA.BeamWidthYError(), spotB.BeamWidthYError(), "witdhY"),
0281        calcSignificance(spotA.sigmaZ(), spotB.sigmaZ(), spotA.sigmaZ0Error(), spotB.sigmaZ0Error(), "witdthZ")}};
0282
0283   if (verbose) {
0284     edm::LogPrint("BeamSpotCompatibilityChecker") << std::string(102, '-');
0285   }
0286
0287   return ret;
0288 }
0289
0290 double BeamSpotCompatibilityChecker::computeBeamSpotCompatibility(const reco::BeamSpot& beamSpot1,
0291                                                                   const reco::BeamSpot& beamSpot2) {
0292   reco::Vertex vertex1(
0293       reco::Vertex::Point(beamSpot1.x0(), beamSpot1.y0(), beamSpot1.z0()), beamSpot1.rotatedCovariance3D(), 0, 0, 0);
0294   reco::Vertex vertex2(
0295       reco::Vertex::Point(beamSpot2.x0(), beamSpot2.y0(), beamSpot2.z0()), beamSpot2.rotatedCovariance3D(), 0, 0, 0);
0296
0297   // Calculate distance and significance using VertexDistance3D
0298   VertexDistance3D distanceCalculator;
0299   // double distance = distanceCalculator.distance(vertex1, vertex2).value();  // Euclidean distance
0300   double significance = distanceCalculator.distance(vertex1, vertex2).significance();  // Distance significance
0301
0302   // Return the significance as a measure of compatibility
0303   return significance;  // Smaller values indicate higher compatibility
0304 }
0305
0306 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0307 void BeamSpotCompatibilityChecker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0308   edm::ParameterSetDescription desc;
0309   desc.add<double>("warningThr", 1.)->setComment("Threshold on the signficances to emit a warning");
0310   desc.add<double>("errorThr", 3.)->setComment("Threshold on the signficances to abort the job");
0311   desc.addUntracked<bool>("verbose", false)->setComment("verbose output");
0312   desc.add<edm::InputTag>("bsFromFile", edm::InputTag(""))
0313       ->setComment("edm::InputTag on the BeamSpot from the File (Reference)");
0314
0315   // Conditional parameters based on dbFromEvent
0316   desc.ifValue(
0317       edm::ParameterDescription<bool>(
0318           "dbFromEvent",
0319           true,
0320           true,
0321           edm::Comment("Switch to take the (target) DB beamspot from the event instead of the EventSetup")),
0322       true >> edm::ParameterDescription<edm::InputTag>(
0323                   "bsFromDB",
0324                   edm::InputTag(""),
0325                   true,
0326                   edm::Comment("edm::InputTag on the BeamSpot from the Event (Target) (used if dbFromEvent = true")) or
0327           false >> edm::ParameterDescription<bool>(
0328                        "useTransientRecord",
0329                        false,
0330                        true,
0331                        edm::Comment("Use transient BeamSpot record (used if dbFromEvent = false)")));
0332
0333   descriptions.addWithDefaultLabel(desc);
0334 }
0335
0336 //define this as a plug-in
0337 DEFINE_FWK_MODULE(BeamSpotCompatibilityChecker);