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 event
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 "DataFormats/BeamSpot/interface/BeamSpot.h"
0026 #include "DataFormats/Common/interface/Handle.h"
0027 #include "DataFormats/Scalers/interface/BeamSpotOnline.h"
0028 #include "FWCore/Framework/interface/ESHandle.h"
0029 #include "FWCore/Framework/interface/ESWatcher.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/EventSetup.h"
0032 #include "FWCore/Framework/interface/Frameworkfwd.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/Framework/interface/global/EDAnalyzer.h"
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0037 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0039 #include "FWCore/Utilities/interface/ESGetToken.h"
0040 #include "FWCore/Utilities/interface/InputTag.h"
0041 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0042 //
0043 // class declaration
0044 //
0045
0046 // ancillary struct to check compatibility
0047 namespace {
0048
0049   struct Significance {
0050     Significance(const double& a, const double& b, const double& errA, const double& errB, const std::string& var)
0051         : m_A(a), m_B(b), m_ErrA(errA), m_ErrB(errB), m_var(var) {
0052       if (m_ErrA == 0 && m_ErrB == 0) {
0053         edm::LogError("LogicalError") << "Can't calculate significance when both errors are zero!" << std::endl;
0054       }
0055       m_combinedError = std::sqrt((m_ErrA * m_ErrA) + (m_ErrB * m_ErrB));
0056     }
0057
0058     float getSig(const bool verbose) {
0059       if (verbose) {
0060         edm::LogPrint("BeamSpotCompatibilityChecker")
0061             << std::fixed << std::setprecision(6)  // Set fixed-point format with 3 decimal places
0062             << m_var << ": A= " << std::setw(10) << m_A << " +/- " << std::setw(10) << m_ErrA
0063             << "    B= " << std::setw(5) << m_B << " +/- " << std::setw(10) << m_ErrB
0064             << "    | delta= " << std::setw(10) << std::abs(m_A - m_B) << " +/- " << std::setw(10) << m_combinedError
0065             << "    Sig= " << std::setw(10) << std::abs(m_A - m_B) / m_combinedError << std::endl;
0066       }
0067       return std::abs(m_A - m_B) / m_combinedError;
0068     }
0069
0070   private:
0071     double m_A;
0072     double m_B;
0073     double m_ErrA;
0074     double m_ErrB;
0075     double m_ErrAB;
0076     std::string m_var;
0077     double m_combinedError;
0078   };
0079 }  // namespace
0080
0081 class BeamSpotCompatibilityChecker : public edm::global::EDAnalyzer<> {
0082 public:
0083   explicit BeamSpotCompatibilityChecker(const edm::ParameterSet&);
0084   ~BeamSpotCompatibilityChecker() override = default;
0085
0086   void analyze(edm::StreamID, edm::Event const&, edm::EventSetup const&) const final;
0087   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0088   static std::array<float, 6> compareBS(const reco::BeamSpot& BSA, const reco::BeamSpot& BSB, const bool verbose);
0089   static double computeBeamSpotCompatibility(const reco::BeamSpot& beamSpot1, const reco::BeamSpot& beamSpot2);
0090
0091 private:
0092   // ----------member data ---------------------------
0093   static constexpr int cmToum = 10000;
0094   const double warningThreshold_;
0095   const double throwingThreshold_;
0096   const bool verbose_;
0097   const bool dbFromEvent_;
0098   const edm::EDGetTokenT<reco::BeamSpot> bsFromEventToken_;             // beamSpot from the event
0099   edm::ESGetToken<BeamSpotObjects, BeamSpotObjectsRcd> bsFromDBToken_;  // beamSpot from the DB
0100   edm::EDGetTokenT<reco::BeamSpot> dbBSfromEventToken_;                 // beamSpot from the DB (via event)
0101 };
0102
0103 //
0104 // constructors and destructor
0105 //
0106 BeamSpotCompatibilityChecker::BeamSpotCompatibilityChecker(const edm::ParameterSet& iConfig)
0107     : warningThreshold_(iConfig.getParameter<double>("warningThr")),
0108       throwingThreshold_(iConfig.getParameter<double>("errorThr")),
0109       verbose_(iConfig.getUntrackedParameter<bool>("verbose", false)),
0110       dbFromEvent_(iConfig.getParameter<bool>("dbFromEvent")),
0111       bsFromEventToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("bsFromEvent"))) {
0112   //now do what ever initialization is needed
0113   if (warningThreshold_ > throwingThreshold_) {
0114     throw cms::Exception("ConfigurationError")
0115         << __PRETTY_FUNCTION__ << "\n Warning threshold (" << warningThreshold_
0116         << ") cannot be smaller than the throwing threshold (" << throwingThreshold_ << ")" << std::endl;
0117   }
0118   if (dbFromEvent_) {
0119     edm::LogWarning("BeamSpotCompatibilityChecker")
0120         << "!!!! Warning !!!\nThe Database Beam Spot is going to be taken from the Event via BeamSpotProducer!";
0121     dbBSfromEventToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("bsFromDB"));
0122   } else {
0123     bsFromDBToken_ = esConsumes<BeamSpotObjects, BeamSpotObjectsRcd>();
0124   }
0125 }
0126
0127 //
0128 // member functions
0129 //
0130
0131 // ------------ method called for each event  ------------
0132 void BeamSpotCompatibilityChecker::analyze(edm::StreamID sid,
0133                                            edm::Event const& iEvent,
0134                                            edm::EventSetup const& iSetup) const {
0135   using namespace edm;
0136   reco::BeamSpot spotEvent, spotDB;
0137
0138   edm::Handle<reco::BeamSpot> beamSpotFromEventHandle;
0139   iEvent.getByToken(bsFromEventToken_, beamSpotFromEventHandle);
0140   spotEvent = *beamSpotFromEventHandle;
0141
0142   double evt_BSx0 = spotEvent.x0();
0143   double evt_BSy0 = spotEvent.y0();
0144   double evt_BSz0 = spotEvent.z0();
0145   double evt_Beamsigmaz = spotEvent.sigmaZ();
0146   double evt_BeamWidthX = spotEvent.BeamWidthX();
0147   double evt_BeamWidthY = spotEvent.BeamWidthY();
0148
0149   if (!dbFromEvent_) {
0150     edm::ESHandle<BeamSpotObjects> beamhandle = iSetup.getHandle(bsFromDBToken_);
0151     const BeamSpotObjects* aSpot = beamhandle.product();
0152
0153     // translate from BeamSpotObjects to reco::BeamSpot
0154     reco::BeamSpot::Point apoint(aSpot->x(), aSpot->y(), aSpot->z());
0155
0156     reco::BeamSpot::CovarianceMatrix matrix;
0157     for (int i = 0; i < reco::BeamSpot::dimension; ++i) {
0158       for (int j = 0; j < reco::BeamSpot::dimension; ++j) {
0159         matrix(i, j) = aSpot->covariance(i, j);
0160       }
0161     }
0162
0163     // this assume beam width same in x and y
0164     spotDB = reco::BeamSpot(apoint, aSpot->sigmaZ(), aSpot->dxdz(), aSpot->dydz(), aSpot->beamWidthX(), matrix);
0165   } else {
0166     // take the DB beamspot from the event (different label)
0167     edm::Handle<reco::BeamSpot> beamSpotFromDBHandle;
0168     iEvent.getByToken(dbBSfromEventToken_, beamSpotFromDBHandle);
0169     spotDB = *beamSpotFromDBHandle;
0170   }
0171
0172   double db_BSx0 = spotDB.x0();
0173   double db_BSy0 = spotDB.y0();
0174   double db_BSz0 = spotDB.z0();
0175   double db_Beamsigmaz = spotDB.sigmaZ();
0176   double db_BeamWidthX = spotDB.BeamWidthX();
0177   double db_BeamWidthY = spotDB.BeamWidthY();
0178
0179   double deltaX0 = evt_BSx0 - db_BSx0;
0180   double deltaY0 = evt_BSy0 - db_BSy0;
0181   double deltaZ0 = evt_BSz0 - db_BSz0;
0182   double deltaSigmaX = evt_BeamWidthX - db_BeamWidthX;
0183   double deltaSigmaY = evt_BeamWidthY - db_BeamWidthY;
0184   double deltaSigmaZ = evt_Beamsigmaz - db_Beamsigmaz;
0185
0186   if (verbose_) {
0187     edm::LogPrint("BeamSpotCompatibilityChecker") << "BS from Event: \n" << spotEvent << std::endl;
0188     edm::LogPrint("BeamSpotCompatibilityChecker") << "BS from DB: \n" << spotDB << std::endl;
0189   }
0190
0191   auto significances = compareBS(spotDB, spotEvent, verbose_);
0192   std::vector<std::string> labels = {"x0", "y0", "z0", "sigmaX", "sigmaY", "sigmaZ"};
0193
0194   std::string msg = " |delta X0|=" + std::to_string(std::abs(deltaX0) * cmToum) +
0195                     " um |delta Y0|=" + std::to_string(std::abs(deltaY0) * cmToum) +
0196                     " um |delta Z0|=" + std::to_string(std::abs(deltaZ0) * cmToum) +
0197                     " um |delta sigmaX|=" + std::to_string(std::abs(deltaSigmaX) * cmToum) +
0198                     " um |delta sigmaY|=" + std::to_string(std::abs(deltaSigmaY) * cmToum) +
0199                     " um |delta sigmaZ|=" + std::to_string(std::abs(deltaSigmaZ)) + " cm";
0200   if (verbose_) {
0201     edm::LogPrint("BeamSpotCompatibilityChecker") << msg.c_str() << std::endl;
0202   }
0203
0204   for (unsigned int i = 0; i < 3; i++) {
0205     auto sig = significances.at(i);
0206     if (sig > throwingThreshold_) {
0207       edm::LogError("BeamSpotCompatibilityChecker") << msg.c_str() << std::endl;
0208       throw cms::Exception("BeamSpotCompatibilityChecker")
0209           << "[" << __PRETTY_FUNCTION__ << "] \n DB-Event BeamSpot " << labels.at(i) << " distance sigificance " << sig
0210           << ", exceeds the threshold of " << throwingThreshold_ << "!" << std::endl;
0211     } else if (sig > warningThreshold_) {
0212       edm::LogWarning("BeamSpotCompatibilityChecker") << msg.c_str() << std::endl;
0213       edm::LogWarning("BeamSpotCompatibilityChecker")
0214           << "[" << __PRETTY_FUNCTION__ << "]  \n  DB-Event BeamSpot " << labels.at(i) << " distance significance "
0215           << sig << ", exceeds the threshold of " << warningThreshold_ << "!" << std::endl;
0216     }
0217   }
0218 }
0219
0220 std::array<float, 6> BeamSpotCompatibilityChecker::compareBS(const reco::BeamSpot& spotA,
0221                                                              const reco::BeamSpot& spotB,
0222                                                              const bool verbose) {
0223   // Lambda to calculate the significance
0224   auto calcSignificance = [&](auto a, auto b, auto aErr, auto bErr, auto var) {
0225     return Significance(a, b, aErr, bErr, var).getSig(verbose);
0226   };
0227
0228   // Populate the array using the lambda
0229   std::array<float, 6> ret = {
0230       {calcSignificance(spotA.x0(), spotB.x0(), spotA.x0Error(), spotB.x0Error(), "x"),
0231        calcSignificance(spotA.y0(), spotB.y0(), spotA.y0Error(), spotB.y0Error(), "y"),
0232        calcSignificance(spotA.z0(), spotB.z0(), spotA.z0Error(), spotB.z0Error(), "z"),
0233        calcSignificance(
0234            spotA.BeamWidthX(), spotB.BeamWidthX(), spotA.BeamWidthXError(), spotB.BeamWidthXError(), "widthX"),
0235        calcSignificance(
0236            spotA.BeamWidthY(), spotB.BeamWidthY(), spotA.BeamWidthYError(), spotB.BeamWidthYError(), "witdhY"),
0237        calcSignificance(spotA.sigmaZ(), spotB.sigmaZ(), spotA.sigmaZ0Error(), spotB.sigmaZ0Error(), "witdthZ")}};
0238
0239   return ret;
0240 }
0241
0242 double BeamSpotCompatibilityChecker::computeBeamSpotCompatibility(const reco::BeamSpot& beamSpot1,
0243                                                                   const reco::BeamSpot& beamSpot2) {
0244   reco::Vertex vertex1(
0245       reco::Vertex::Point(beamSpot1.x0(), beamSpot1.y0(), beamSpot1.z0()), beamSpot1.rotatedCovariance3D(), 0, 0, 0);
0246   reco::Vertex vertex2(
0247       reco::Vertex::Point(beamSpot2.x0(), beamSpot2.y0(), beamSpot2.z0()), beamSpot2.rotatedCovariance3D(), 0, 0, 0);
0248
0249   // Calculate distance and significance using VertexDistance3D
0250   VertexDistance3D distanceCalculator;
0251   // double distance = distanceCalculator.distance(vertex1, vertex2).value();  // Euclidean distance
0252   double significance = distanceCalculator.distance(vertex1, vertex2).significance();  // Distance significance
0253
0254   // Return the significance as a measure of compatibility
0255   return significance;  // Smaller values indicate higher compatibility
0256 }
0257
0258 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0259 void BeamSpotCompatibilityChecker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0260   edm::ParameterSetDescription desc;
0261   desc.add<double>("warningThr", 1.)->setComment("Threshold on the signficances to emit a warning");
0262   desc.add<double>("errorThr", 3.)->setComment("Threshold on the signficances to abort the job");
0263   desc.addUntracked<bool>("verbose", false)->setComment("verbose output");
0264   desc.add<edm::InputTag>("bsFromEvent", edm::InputTag(""))
0265       ->setComment("edm::InputTag on the BeamSpot from the Event (Reference)");
0266   desc.add<bool>("dbFromEvent", false)
0267       ->setComment("Switch to take the (target) DB beamspot from the event instead of the EventSetup");
0268   desc.add<edm::InputTag>("bsFromDB", edm::InputTag(""))
0269       ->setComment("edm::InputTag on the BeamSpot from the Event (Target)\n To be used only if dbFromEvent is True");
0270   descriptions.addWithDefaultLabel(desc);
0271 }
0272
0273 //define this as a plug-in
0274 DEFINE_FWK_MODULE(BeamSpotCompatibilityChecker);