Merge pull request #45637 from mmusich/mm_add_BeamSpotCompatibilityCh…

…ecker_14_0_X

[14.0.X] Introduce `BeamSpotCompatibilityChecker` and corresponding tests

RecoVertex/BeamSpotProducer/plugins/BeamSpotCompatibilityChecker.cc

@@ -0,0 +1,248 @@
+// -*- C++ -*-
+//
+// Package:    RecoVertex/BeamSpotProducer
+// Class:      BeamSpotCompatibilityChecker
+//
+/**\class BeamSpotCompatibilityChecker BeamSpotCompatibilityChecker.cc RecoVertex/BeamSpotProducer/plugins/BeamSpotCompatibilityChecker.cc
+
+ Description: Class to check the compatibilty between the BeamSpot payload in the database and the one in the event
+
+ Implementation:
+     Makes use of the Significance struct to establish how compatible are the data members of the two BeamSpots in input
+*/
+//
+// Original Author:  Marco Musich
+//         Created:  Thu, 23 Apr 2020 09:00:45 GMT
+//
+//
+
+// system include files
+#include <memory>
+
+// user include files
+#include "CondFormats/BeamSpotObjects/interface/BeamSpotObjects.h"
+#include "CondFormats/DataRecord/interface/BeamSpotObjectsRcd.h"
+#include "DataFormats/BeamSpot/interface/BeamSpot.h"
+#include "DataFormats/Common/interface/Handle.h"
+#include "DataFormats/Scalers/interface/BeamSpotOnline.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/ESWatcher.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/global/EDAnalyzer.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "FWCore/Utilities/interface/ESGetToken.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+
+//
+// class declaration
+//
+
+// ancillary struct to check compatibility
+namespace {
+
+  struct Significance {
+    Significance(const double& a, const double& b, const double& errA, const double& errB)
+        : m_A(a), m_B(b), m_ErrA(errA), m_ErrB(errB) {
+      if (m_ErrA == 0 && m_ErrB == 0) {
+        edm::LogError("LogicalError") << "Can't calculate significance when both errors are zero!" << std::endl;
+      }
+      m_combinedError = std::sqrt((m_ErrA * m_ErrA) + (m_ErrB * m_ErrB));
+    }
+
+    float getSig(bool verbose) {
+      if (verbose) {
+        edm::LogInfo("BeamSpotCompatibilityChecker")
+            << "A= " << m_A << "+/-" << m_ErrA << "    "
+            << "B= " << m_B << "+/-" << m_ErrB << " | delta=" << std::abs(m_A - m_B) << "+/-" << m_combinedError
+            << std::endl;
+      }
+      return std::abs(m_A - m_B) / m_combinedError;
+    }
+
+  private:
+    double m_A;
+    double m_B;
+    double m_ErrA;
+    double m_ErrB;
+    double m_ErrAB;
+    double m_combinedError;
+  };
+}  // namespace
+
+class BeamSpotCompatibilityChecker : public edm::global::EDAnalyzer<> {
+public:
+  explicit BeamSpotCompatibilityChecker(const edm::ParameterSet&);
+  ~BeamSpotCompatibilityChecker() override = default;
+
+  void analyze(edm::StreamID, edm::Event const&, edm::EventSetup const&) const final;
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+  static std::array<float, 6> compareBS(const reco::BeamSpot& BSA, const reco::BeamSpot& BSB);
+
+private:
+  // ----------member data ---------------------------
+  static constexpr int cmToum = 10000;
+  const double warningThreshold_;
+  const double throwingThreshold_;
+  const bool verbose_;
+  const bool dbFromEvent_;
+  const edm::EDGetTokenT<reco::BeamSpot> bsFromEventToken_;             // beamSpot from the event
+  edm::ESGetToken<BeamSpotObjects, BeamSpotObjectsRcd> bsFromDBToken_;  // beamSpot from the DB
+  edm::EDGetTokenT<reco::BeamSpot> dbBSfromEventToken_;                 // beamSpot from the DB (via event)
+};
+
+//
+// constructors and destructor
+//
+BeamSpotCompatibilityChecker::BeamSpotCompatibilityChecker(const edm::ParameterSet& iConfig)
+    : warningThreshold_(iConfig.getParameter<double>("warningThr")),
+      throwingThreshold_(iConfig.getParameter<double>("errorThr")),
+      verbose_(iConfig.getUntrackedParameter<bool>("verbose", false)),
+      dbFromEvent_(iConfig.getParameter<bool>("dbFromEvent")),
+      bsFromEventToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("bsFromEvent"))) {
+  //now do what ever initialization is needed
+  if (warningThreshold_ > throwingThreshold_) {
+    throw cms::Exception("ConfigurationError")
+        << __PRETTY_FUNCTION__ << "\n Warning threshold (" << warningThreshold_
+        << ") cannot be smaller than the throwing threshold (" << throwingThreshold_ << ")" << std::endl;
+  }
+  if (dbFromEvent_) {
+    edm::LogWarning("BeamSpotCompatibilityChecker")
+        << "!!!! Warning !!!\nThe Database Beam Spot is going to be taken from the Event via BeamSpotProducer!";
+    dbBSfromEventToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("bsFromDB"));
+  } else {
+    bsFromDBToken_ = esConsumes<BeamSpotObjects, BeamSpotObjectsRcd>();
+  }
+}
+
+//
+// member functions
+//
+
+// ------------ method called for each event  ------------
+void BeamSpotCompatibilityChecker::analyze(edm::StreamID sid,
+                                           edm::Event const& iEvent,
+                                           edm::EventSetup const& iSetup) const {
+  using namespace edm;
+  reco::BeamSpot spotEvent, spotDB;
+
+  edm::Handle<reco::BeamSpot> beamSpotFromEventHandle;
+  iEvent.getByToken(bsFromEventToken_, beamSpotFromEventHandle);
+  spotEvent = *beamSpotFromEventHandle;
+
+  double evt_BSx0 = spotEvent.x0();
+  double evt_BSy0 = spotEvent.y0();
+  double evt_BSz0 = spotEvent.z0();
+  double evt_Beamsigmaz = spotEvent.sigmaZ();
+  double evt_BeamWidthX = spotEvent.BeamWidthX();
+  double evt_BeamWidthY = spotEvent.BeamWidthY();
+
+  if (!dbFromEvent_) {
+    edm::ESHandle<BeamSpotObjects> beamhandle = iSetup.getHandle(bsFromDBToken_);
+    const BeamSpotObjects* aSpot = beamhandle.product();
+
+    // translate from BeamSpotObjects to reco::BeamSpot
+    reco::BeamSpot::Point apoint(aSpot->x(), aSpot->y(), aSpot->z());
+
+    reco::BeamSpot::CovarianceMatrix matrix;
+    for (int i = 0; i < reco::BeamSpot::dimension; ++i) {
+      for (int j = 0; j < reco::BeamSpot::dimension; ++j) {
+        matrix(i, j) = aSpot->covariance(i, j);
+      }
+    }
+
+    // this assume beam width same in x and y
+    spotDB = reco::BeamSpot(apoint, aSpot->sigmaZ(), aSpot->dxdz(), aSpot->dydz(), aSpot->beamWidthX(), matrix);
+  } else {
+    // take the DB beamspot from the event (different label)
+    edm::Handle<reco::BeamSpot> beamSpotFromDBHandle;
+    iEvent.getByToken(dbBSfromEventToken_, beamSpotFromDBHandle);
+    spotDB = *beamSpotFromDBHandle;
+  }
+
+  double db_BSx0 = spotDB.x0();
+  double db_BSy0 = spotDB.y0();
+  double db_BSz0 = spotDB.z0();
+  double db_Beamsigmaz = spotDB.sigmaZ();
+  double db_BeamWidthX = spotDB.BeamWidthX();
+  double db_BeamWidthY = spotDB.BeamWidthY();
+
+  double deltaX0 = evt_BSx0 - db_BSx0;
+  double deltaY0 = evt_BSy0 - db_BSy0;
+  double deltaZ0 = evt_BSz0 - db_BSz0;
+  double deltaSigmaX = evt_BeamWidthX - db_BeamWidthX;
+  double deltaSigmaY = evt_BeamWidthY - db_BeamWidthY;
+  double deltaSigmaZ = evt_Beamsigmaz - db_Beamsigmaz;
+
+  if (verbose_) {
+    edm::LogPrint("BeamSpotCompatibilityChecker") << "BS from Event: \n" << spotEvent << std::endl;
+    edm::LogPrint("BeamSpotCompatibilityChecker") << "BS from DB: \n" << spotDB << std::endl;
+  }
+
+  auto significances = compareBS(spotDB, spotEvent);
+  std::vector<std::string> labels = {"x0", "y0", "z0", "sigmaX", "sigmaY", "sigmaZ"};
+
+  std::string msg = " |delta X0|=" + std::to_string(std::abs(deltaX0) * cmToum) +
+                    " um |delta Y0|=" + std::to_string(std::abs(deltaY0) * cmToum) +
+                    " um |delta Z0|=" + std::to_string(std::abs(deltaZ0) * cmToum) +
+                    " um |delta sigmaX|=" + std::to_string(std::abs(deltaSigmaX) * cmToum) +
+                    " um |delta sigmaY|=" + std::to_string(std::abs(deltaSigmaY) * cmToum) +
+                    " um |delta sigmaZ|=" + std::to_string(std::abs(deltaSigmaZ)) + "cm";
+  if (verbose_) {
+    edm::LogPrint("BeamSpotCompatibilityChecker") << msg.c_str() << std::endl;
+  }
+
+  for (unsigned int i = 0; i < 3; i++) {
+    auto sig = significances.at(i);
+    if (sig > throwingThreshold_) {
+      edm::LogError("BeamSpotCompatibilityChecker") << msg.c_str() << std::endl;
+      throw cms::Exception("BeamSpotCompatibilityChecker")
+          << "[" << __PRETTY_FUNCTION__ << "] \n DB-Event BeamSpot " << labels.at(i) << " distance sigificance " << sig
+          << ", exceeds the threshold of " << throwingThreshold_ << "!" << std::endl;
+    } else if (sig > warningThreshold_) {
+      edm::LogWarning("BeamSpotCompatibilityChecker") << msg.c_str() << std::endl;
+      edm::LogWarning("BeamSpotCompatibilityChecker")
+          << "[" << __PRETTY_FUNCTION__ << "]  \n  DB-Event BeamSpot " << labels.at(i) << " distance significance "
+          << sig << ", exceeds the threshold of " << warningThreshold_ << "!" << std::endl;
+    }
+  }
+}
+
+std::array<float, 6> BeamSpotCompatibilityChecker::compareBS(const reco::BeamSpot& beamSpotA,
+                                                             const reco::BeamSpot& beamSpotB) {
+  auto xS = Significance(beamSpotA.x0(), beamSpotB.x0(), beamSpotA.x0Error(), beamSpotB.x0Error());
+  auto yS = Significance(beamSpotA.y0(), beamSpotB.y0(), beamSpotA.y0Error(), beamSpotB.y0Error());
+  auto zS = Significance(beamSpotA.z0(), beamSpotB.z0(), beamSpotA.z0Error(), beamSpotB.z0Error());
+  auto xWS = Significance(
+      beamSpotA.BeamWidthX(), beamSpotB.BeamWidthX(), beamSpotA.BeamWidthXError(), beamSpotB.BeamWidthXError());
+  auto yWS = Significance(
+      beamSpotA.BeamWidthY(), beamSpotB.BeamWidthY(), beamSpotA.BeamWidthYError(), beamSpotB.BeamWidthYError());
+  auto zWS = Significance(beamSpotA.sigmaZ(), beamSpotB.sigmaZ(), beamSpotA.sigmaZ0Error(), beamSpotB.sigmaZ0Error());
+
+  std::array<float, 6> ret = {
+      {xS.getSig(false), yS.getSig(false), zS.getSig(false), xWS.getSig(false), yWS.getSig(false), zWS.getSig(false)}};
+  return ret;
+}
+
+// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
+void BeamSpotCompatibilityChecker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<double>("warningThr", 1.)->setComment("Threshold on the signficances to emit a warning");
+  desc.add<double>("errorThr", 3.)->setComment("Threshold on the signficances to abort the job");
+  desc.addUntracked<bool>("verbose", false)->setComment("verbose output");
+  desc.add<edm::InputTag>("bsFromEvent", edm::InputTag(""))
+      ->setComment("edm::InputTag on the BeamSpot from the Event (Reference)");
+  desc.add<bool>("dbFromEvent", false)
+      ->setComment("Switch to take the (target) DB beamspot from the event instead of the EventSetup");
+  desc.add<edm::InputTag>("bsFromDB", edm::InputTag(""))
+      ->setComment("edm::InputTag on the BeamSpot from the Event (Target)\n To be used only if dbFromEvent is True");
+  descriptions.addWithDefaultLabel(desc);
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(BeamSpotCompatibilityChecker);

RecoVertex/BeamSpotProducer/plugins/BuildFile.xml

@@ -48,6 +48,10 @@
   <flags EDM_PLUGIN="1"/>
 </library>
 
+<library file="BeamSpotCompatibilityChecker.cc" name="BeamSpotCompatibilityChecker">
+  <flags EDM_PLUGIN="1"/>
+</library>
+
 <iftool name="cuda">
   <library file="BeamSpotToCUDA.cc" name="BeamSpotToCUDA">
     <use name="cuda"/>

RecoVertex/BeamSpotProducer/test/BuildFile.xml

@@ -11,7 +11,16 @@
 
 <test name="testReadWriteBSFromDB" command="testReadWriteBSFromDB.sh"/>
 
+<test name="testBeamSpotCompatibilityEventData" command="testBeamSpotCompatibility.sh"/>
+
 <bin file="testBeamSpotAnalyzer.cc">
   <use name="FWCore/TestProcessor"/>
   <use name="catch2"/>
 </bin>
+
+<bin file="testBeamSpotCompatibility.cc">
+  <use name="CondFormats/BeamSpotObjects"/>
+  <use name="CondFormats/DataRecord"/>
+  <use name="FWCore/TestProcessor"/>
+  <use name="catch2"/>
+</bin>

RecoVertex/BeamSpotProducer/test/testBeamSpotCompatibility.cc

@@ -0,0 +1,63 @@
+#include "FWCore/TestProcessor/interface/TestProcessor.h"
+#include "DataFormats/BeamSpot/interface/BeamSpot.h"
+#include "RecoVertex/BeamSpotProducer/plugins/BeamSpotCompatibilityChecker.cc"
+
+#define CATCH_CONFIG_MAIN
+#include "catch.hpp"
+
+//_________________________________________________________
+TEST_CASE("Significance Calculation", "[Significance]") {
+  double a = 10.0;
+  double b = 12.0;
+  double errA = 1.0;
+  double errB = 1.5;
+
+  Significance sig(a, b, errA, errB);
+  float significance = sig.getSig(false);
+
+  // Correct the expected value based on actual calculation
+  REQUIRE(significance == Approx(1.1094).epsilon(10e-6));
+}
+
+//_________________________________________________________
+TEST_CASE("BeamSpot Compatibility Checker", "[compareBS]") {
+  reco::BeamSpot::Point pointA(1.0, 2.0, 3.0);
+  reco::BeamSpot::Point pointB(1.1, 2.1, 3.1);
+
+  reco::BeamSpot::CovarianceMatrix matrixA;
+  reco::BeamSpot::CovarianceMatrix matrixB;
+
+  // Initialize matrices with non-zero values
+  for (int i = 0; i < reco::BeamSpot::dimension; ++i) {
+    for (int j = 0; j < reco::BeamSpot::dimension; ++j) {
+      matrixA(i, j) = 0.01 * (i + 1) * (j + 1);
+      matrixB(i, j) = 0.02 * (i + 1) * (j + 1);
+    }
+  }
+
+  reco::BeamSpot beamSpotA(pointA, 4.0, 0.01, 0.01, 0.1, matrixA);
+  reco::BeamSpot beamSpotB(pointB, 4.2, 0.01, 0.01, 0.12, matrixB);
+
+  // Create the edm::ParameterSet with the required parameters
+  edm::ParameterSet pset;
+  pset.addParameter<double>("warningThr", 1.0);
+  pset.addParameter<double>("errorThr", 2.0);
+  pset.addParameter<edm::InputTag>("bsFromEvent", edm::InputTag(""));
+  pset.addParameter<bool>("dbFromEvent", false);
+  pset.addParameter<edm::InputTag>("bsFromDB", edm::InputTag(""));
+
+  BeamSpotCompatibilityChecker checker(pset);
+  auto significances = checker.compareBS(beamSpotA, beamSpotB);
+
+  // Print significances
+  for (size_t i = 0; i < significances.size(); ++i) {
+    std::cout << "Significance[" << i << "]: " << significances[i] << std::endl;
+  }
+
+  REQUIRE(significances[0] == Approx(0.57735).epsilon(10e-6));    // x0 significance
+  REQUIRE(significances[1] == Approx(0.288675).epsilon(10e-6));   // y0 significance
+  REQUIRE(significances[2] == Approx(0.19245).epsilon(10e-6));    // z0 significance
+  REQUIRE(significances[3] == Approx(0.0164957).epsilon(10e-6));  // sigmaX significance
+  REQUIRE(significances[4] == Approx(0.0164957).epsilon(10e-6));  // sigmaY significance
+  REQUIRE(significances[5] == Approx(0.288675).epsilon(10e-6));   // sigmaZ significance
+}

RecoVertex/BeamSpotProducer/test/testBeamSpotCompatibility.sh

@@ -0,0 +1,33 @@
+#!/bin/sh
+
+check_for_success() {
+    "${@}" && echo -e "\n ---> Passed test of '${@}'\n\n" || exit 1
+}
+
+check_for_failure() {
+    "${@}" && exit 1 || echo -e "\n ---> Passed test of '${@}'\n\n"
+}
+
+echo "TESTING BeamSpot compatibility check code ..."
+
+
+echo "========================================"
+echo "testing failing comparison"
+echo "----------------------------------------"
+echo
+
+check_for_failure cmsRun ${SCRAM_TEST_PATH}/testBeamSpotCompatibility_cfg.py
+
+echo "========================================"
+echo "testing successful comparison (with warning)"
+echo "----------------------------------------"
+echo
+
+check_for_success cmsRun ${SCRAM_TEST_PATH}/testBeamSpotCompatibility_cfg.py errorThreshold=100000
+
+echo "========================================"
+echo "testing failed comparison (DB from event)"
+echo "----------------------------------------"
+echo
+
+check_for_failure cmsRun ${SCRAM_TEST_PATH}/testBeamSpotCompatibility_cfg.py dbFromEvent=True