Update SimBeamSpotObjects to handle different EvtVtxGenerators

CondCore/BeamSpotPlugins/interface/BeamSpotPayloadInspectorHelper.h

@@ -715,17 +715,21 @@ namespace beamSpotPI {
 namespace simBeamSpotPI {
 
   enum parameters {
-    X = 0,              // 0  - Positions
-    Y = 1,              // 1
-    Z = 2,              // 2
-    sigmaZ = 3,         // 3  - Widths
-    betaStar = 4,       // 4
-    emittance = 5,      // 5
-    expTransWidth = 6,  // 6  - from LPC-like calculation
-    phi = 7,            // 7  - Additional parameters
-    alpha = 8,          // 8
-    timeOffset = 9,     // 9
-    END_OF_TYPES = 10,
+    X = 0,            // 0  - Positions
+    Y = 1,            // 1
+    Z = 2,            // 2
+    meanX = 3,        // 3
+    meanY = 4,        // 4
+    meanZ = 5,        // 5
+    sigmaX = 6,       // 6  - Widths
+    sigmaY = 7,       // 7
+    sigmaZ = 8,       // 8
+    betaStar = 9,     // 9
+    emittance = 10,   // 10
+    phi = 11,         // 11  - Additional parameters
+    alpha = 12,       // 12
+    timeOffset = 13,  // 13
+    END_OF_TYPES = 14,
   };
 
   /************************************************/
@@ -737,14 +741,22 @@ namespace simBeamSpotPI {
         return (addUnits ? "Y [cm]" : "Y");
       case Z:
         return (addUnits ? "Z [cm]" : "Z");
+      case meanX:
+        return (addUnits ? "MeanX [cm]" : "meanX");
+      case meanY:
+        return (addUnits ? "MeanY [cm]" : "meanY");
+      case meanZ:
+        return (addUnits ? "MeanZ [cm]" : "meanZ");
+      case sigmaX:
+        return (addUnits ? "#sigma_{X} [#mum]" : "sigmaX");
+      case sigmaY:
+        return (addUnits ? "#sigma_{Y} [#mum]" : "sigmaY");
       case sigmaZ:
         return (addUnits ? "#sigma_{Z} [cm]" : "sigmaZ");
       case betaStar:
         return (addUnits ? "#beta* [cm]" : "BetaStar");
       case emittance:
         return (addUnits ? "Emittance [cm]" : "Emittance");
-      case expTransWidth:
-        return (addUnits ? "#sigma^{trans}_{xy} [#mum]" : "Exp. trans width");
       case phi:
         return (addUnits ? "Phi [rad]" : "Phi");
       case alpha:
@@ -769,9 +781,12 @@ namespace simBeamSpotPI {
     SimBSParamsHelper(const std::shared_ptr<PayloadType>& bs) {
       // fill in the values
       m_values[parameters::X] = bs->x(), m_values[parameters::Y] = bs->y(), m_values[parameters::Z] = bs->z();
-      m_values[parameters::sigmaZ] = bs->sigmaZ(), m_values[parameters::betaStar] = bs->betaStar(),
-      m_values[parameters::emittance] = bs->emittance();
-      m_values[parameters::expTransWidth] = (1 / std::sqrt(2)) * std::sqrt(bs->emittance() * bs->betaStar()) * 10000.f;
+      m_values[parameters::meanX] = bs->meanX(), m_values[parameters::meanY] = bs->meanY();
+      m_values[parameters::meanZ] = bs->meanZ();
+      m_values[parameters::sigmaX] = bs->sigmaX() * 10000.f;
+      m_values[parameters::sigmaY] = bs->sigmaY() * 10000.f;
+      m_values[parameters::sigmaZ] = bs->sigmaZ();
+      m_values[parameters::betaStar] = bs->betaStar(), m_values[parameters::emittance] = bs->emittance();
       m_values[parameters::phi] = bs->phi(), m_values[parameters::alpha] = bs->alpha(),
       m_values[parameters::timeOffset] = bs->timeOffset();
     }
@@ -842,6 +857,16 @@ namespace simBeamSpotPI {
             return m_payload->y();
           case Z:
             return m_payload->z();
+          case meanX:
+            return m_payload->meanX();
+          case meanY:
+            return m_payload->meanY();
+          case meanZ:
+            return m_payload->meanZ();
+          case sigmaX:
+            return m_payload->sigmaX() * cmToUm;
+          case sigmaY:
+            return m_payload->sigmaY() * cmToUm;
           case sigmaZ:
             return m_payload->sigmaZ();
           case betaStar:
@@ -854,8 +879,6 @@ namespace simBeamSpotPI {
             return m_payload->alpha();
           case timeOffset:
             return m_payload->timeOffset();
-          case expTransWidth:
-            return (1 / std::sqrt(2)) * std::sqrt(m_payload->emittance() * m_payload->betaStar()) * cmToUm;
           case END_OF_TYPES:
             return ret;
           default:

CondCore/BeamSpotPlugins/test/testBeamSpotPayloadInspector.cpp

@@ -84,7 +84,7 @@ int main(int argc, char** argv) {
   // SimBeamSpot
   std::string prepConnectionString("frontier://FrontierPrep/CMS_CONDITIONS");
 
-  tag = "SimBeamSpotObjects_Realistic25ns13p6TeVEOY2022Collision_v0_mc";
+  tag = "SimBeamSpot_Realistic25ns13p6TeVEOY2022CollisionVtxSmearingParameters_v1_mc";
   start = static_cast<unsigned long long>(1);
   end = static_cast<unsigned long long>(1);
 
@@ -98,8 +98,8 @@ int main(int argc, char** argv) {
   histoSimParametersDiff.process(prepConnectionString, PI::mk_input(tag, start, end));
   edm::LogPrint("testBeamSpotPayloadInspector") << histoSimParametersDiff.data() << std::endl;
 
-  tag1 = "SimBeamSpotObjects_Realistic25ns13p6TeVEOY2022Collision_v0_mc";
-  tag2 = "SimBeamSpotObjects_Realistic25ns13p6TeVEarly2023Collision_v0_mc";
+  tag1 = "SimBeamSpot_Realistic25ns13p6TeVEOY2022CollisionVtxSmearingParameters_v1_mc";
+  tag2 = "SimBeamSpot_Realistic25ns13p6TeVEarly2023CollisionVtxSmearingParameters_v1_mc";
   start = static_cast<unsigned long long>(1);
 
   SimBeamSpotParametersDiffTwoTags histoSimParametersDiffTwoTags;

CondFormats/BeamSpotObjects/interface/SimBeamSpotObjects.h

@@ -3,7 +3,21 @@
 
 /** \class SimBeamSpotObjects
  *
- * provide the vertex smearing parameters from DB
+ * Provide the vertex smearing parameters from DB
+ *
+ * This Object contains the parameters needed by the vtx smearing functions used up to Run 3:
+ *   - BetafuncEvtVtxGenerator (realistic Run1/Run2/Run3 conditions)
+ *     Parameters used:
+ *        - fX0, fY0, fZ0
+ *        - fSigmaZ
+ *        - fbetastar, femittance
+ *        - fPhi, fAlpha
+ *        - fTimeOffset
+ *   - GaussEvtVtxGenerator (design Run1/Run2/Run3 conditions)
+ *     Parameters used:
+ *        - fMeanX, fMeanY, fMeanZ
+ *        - fSigmaX, fSigmaY, fSigmaZ
+ *        - fTimeOffset
  *
  */
 
@@ -18,6 +32,11 @@ class SimBeamSpotObjects {
     fX0 = 0.0;
     fY0 = 0.0;
     fZ0 = 0.0;
+    fMeanX = 0.0;
+    fMeanY = 0.0;
+    fMeanZ = 0.0;
+    fSigmaX = -1.0;
+    fSigmaY = -1.0;
     fSigmaZ = 0.0;
     fbetastar = 0.0;
     femittance = 0.0;
@@ -32,7 +51,13 @@ class SimBeamSpotObjects {
   void setX(double val) { fX0 = val; }
   void setY(double val) { fY0 = val; }
   void setZ(double val) { fZ0 = val; }
-  /// set sigmaZ
+  /// set meanX, meanY, meanZ
+  void setMeanX(double val) { fMeanX = val; }
+  void setMeanY(double val) { fMeanY = val; }
+  void setMeanZ(double val) { fMeanZ = val; }
+  /// set sigmaX, sigmaY, sigmaZ
+  void setSigmaX(double val) { fSigmaX = val; }
+  void setSigmaY(double val) { fSigmaY = val; }
   void setSigmaZ(double val) { fSigmaZ = val; }
   /// set BetaStar and Emittance
   void setBetaStar(double val) { fbetastar = val; }
@@ -42,31 +67,33 @@ class SimBeamSpotObjects {
   void setAlpha(double val) { fAlpha = val; }
   void setTimeOffset(double val) { fTimeOffset = val; }
 
-  /// get X position
+  /// get X, Y, Z position
   double x() const { return fX0; }
-  /// get Y position
   double y() const { return fY0; }
-  /// get Z position
   double z() const { return fZ0; }
-  /// get sigmaZ
+  /// get meanX, meanY, meanZ position
+  double meanX() const { return fMeanX; }
+  double meanY() const { return fMeanY; }
+  double meanZ() const { return fMeanZ; }
+  /// get sigmaX, sigmaY, sigmaZ
+  double sigmaX() const;
+  double sigmaY() const;
   double sigmaZ() const { return fSigmaZ; }
-  /// get BetaStar
+  /// get BetaStar and Emittance
   double betaStar() const { return fbetastar; }
-  /// get Emittance
   double emittance() const { return femittance; }
-  /// get Phi
+  /// get Phi, Alpha and TimeOffset
   double phi() const { return fPhi; }
-  /// get Alpha
   double alpha() const { return fAlpha; }
-  /// get TimeOffset
   double timeOffset() const { return fTimeOffset; }
 
   /// print sim beam spot parameters
   void print(std::stringstream& ss) const;
 
 private:
   double fX0, fY0, fZ0;
-  double fSigmaZ;
+  double fMeanX, fMeanY, fMeanZ;
+  double fSigmaX, fSigmaY, fSigmaZ;
   double fbetastar, femittance;
   double fPhi, fAlpha;
   double fTimeOffset;

CondFormats/BeamSpotObjects/src/SimBeamSpotObjects.cc

@@ -1,13 +1,37 @@
 #include "CondFormats/BeamSpotObjects/interface/SimBeamSpotObjects.h"
 
+#include <cmath>
 #include <iostream>
 
+// Get SigmaX and SigmaY:
+// - directly fSigmaX if >= 0 (in case of Gaussian Smearing)
+// - else from LPC-like calculation (in case of BetaFunc Smearing)
+double SimBeamSpotObjects::sigmaX() const {
+  if (fSigmaX >= 0.)  // Gaussian smearing
+    return fSigmaX;
+  else  // BetaFunc smearing
+    return (1 / std::sqrt(2)) * std::sqrt(femittance * fbetastar);
+}
+
+double SimBeamSpotObjects::sigmaY() const {
+  if (fSigmaY >= 0.)  // Gaussian smearing
+    return fSigmaY;
+  else  // BetaFunc smearing
+    return (1 / std::sqrt(2)) * std::sqrt(femittance * fbetastar);
+}
+
+// Printout SimBeamSpotObjects
 void SimBeamSpotObjects::print(std::stringstream& ss) const {
   ss << "-----------------------------------------------------\n"
      << "              Sim Beam Spot Data\n\n"
      << "       X0     = " << x() << " [cm]\n"
      << "       Y0     = " << y() << " [cm]\n"
      << "       Z0     = " << z() << " [cm]\n"
+     << "    MeanX     = " << meanX() << " [cm]\n"
+     << "    MeanY     = " << meanY() << " [cm]\n"
+     << "    MeanZ     = " << meanZ() << " [cm]\n"
+     << " Sigma X0     = " << sigmaX() << " [cm]\n"
+     << " Sigma Y0     = " << sigmaY() << " [cm]\n"
      << " Sigma Z0     = " << sigmaZ() << " [cm]\n"
      << " Beta star    = " << betaStar() << " [cm]\n"
      << " Emittance X  = " << emittance() << " [cm]\n"

CondTools/BeamSpot/plugins/BeamProfile2DBReader.cc

@@ -56,7 +56,8 @@ class BeamProfile2DBReader : public edm::one::EDAnalyzer<edm::one::SharedResourc
     int run;
     int ls;
     double fX0, fY0, fZ0;
-    double fSigmaZ;
+    double fMeanX, fMeanY, fMeanZ;
+    double fSigmaX, fSigmaY, fSigmaZ;
     double fbetastar, femittance;
     double fPhi, fAlpha;
     double fTimeOffset;
@@ -102,6 +103,11 @@ void BeamProfile2DBReader::TheBSfromDB::init() {
   fX0 = dummy_double;
   fY0 = dummy_double;
   fZ0 = dummy_double;
+  fMeanX = dummy_double;
+  fMeanY = dummy_double;
+  fMeanZ = dummy_double;
+  fSigmaX = dummy_double;
+  fSigmaY = dummy_double;
   fSigmaZ = dummy_double;
   fbetastar = dummy_double;
   femittance = dummy_double;
@@ -130,6 +136,11 @@ void BeamProfile2DBReader::analyze(const edm::Event& iEvent, const edm::EventSet
     theBSfromDB_.fX0 = mybeamspot->x();
     theBSfromDB_.fY0 = mybeamspot->y();
     theBSfromDB_.fZ0 = mybeamspot->z();
+    theBSfromDB_.fMeanX = mybeamspot->meanX();
+    theBSfromDB_.fMeanY = mybeamspot->meanY();
+    theBSfromDB_.fMeanZ = mybeamspot->meanZ();
+    theBSfromDB_.fSigmaX = mybeamspot->sigmaX();
+    theBSfromDB_.fSigmaY = mybeamspot->sigmaY();
     theBSfromDB_.fSigmaZ = mybeamspot->sigmaZ();
     theBSfromDB_.fbetastar = mybeamspot->betaStar();
     theBSfromDB_.femittance = mybeamspot->emittance();
@@ -157,6 +168,11 @@ void BeamProfile2DBReader::beginJob() {
   bstree_->Branch("BSx0", &theBSfromDB_.fX0, "BSx0/F");
   bstree_->Branch("BSy0", &theBSfromDB_.fY0, "BSy0/F");
   bstree_->Branch("BSz0", &theBSfromDB_.fZ0, "BSz0/F");
+  bstree_->Branch("BSmeanX", &theBSfromDB_.fMeanX, "BSmeanX/F");
+  bstree_->Branch("BSmeanY", &theBSfromDB_.fMeanY, "BSmeanY/F");
+  bstree_->Branch("BSmeanZ", &theBSfromDB_.fMeanZ, "BSmeanZ/F");
+  bstree_->Branch("Beamsigmax", &theBSfromDB_.fSigmaX, "Beamsigmax/F");
+  bstree_->Branch("Beamsigmay", &theBSfromDB_.fSigmaY, "Beamsigmay/F");
   bstree_->Branch("Beamsigmaz", &theBSfromDB_.fSigmaZ, "Beamsigmaz/F");
   bstree_->Branch("BetaStar", &theBSfromDB_.fbetastar, "BetaStar/F");
   bstree_->Branch("Emittance", &theBSfromDB_.femittance, "Emittance/F");

CondTools/BeamSpot/plugins/BeamProfile2DBWriter.cc

@@ -60,6 +60,11 @@ BeamProfile2DBWriter::BeamProfile2DBWriter(const edm::ParameterSet& iConfig)
   beamSpot_.setX(iConfig.getParameter<double>("X0"));
   beamSpot_.setY(iConfig.getParameter<double>("Y0"));
   beamSpot_.setZ(iConfig.getParameter<double>("Z0"));
+  beamSpot_.setMeanX(iConfig.getParameter<double>("MeanX"));
+  beamSpot_.setMeanY(iConfig.getParameter<double>("MeanY"));
+  beamSpot_.setMeanZ(iConfig.getParameter<double>("MeanZ"));
+  beamSpot_.setSigmaX(iConfig.getParameter<double>("SigmaX"));
+  beamSpot_.setSigmaY(iConfig.getParameter<double>("SigmaY"));
   beamSpot_.setSigmaZ(iConfig.getParameter<double>("SigmaZ"));
   beamSpot_.setAlpha(iConfig.getParameter<double>("Alpha"));
   beamSpot_.setPhi(iConfig.getParameter<double>("Phi"));
@@ -88,12 +93,17 @@ void BeamProfile2DBWriter::fillDescriptions(edm::ConfigurationDescriptions& desc
   edm::ParameterSetDescription desc;
   desc.add<std::string>("recordName", "SimBeamSpotObjectsRcd")
       ->setComment("name of the record to use for the PoolDBOutputService");
-  desc.add<double>("X0")->setComment("in cm");
-  desc.add<double>("Y0")->setComment("in cm");
-  desc.add<double>("Z0")->setComment("in cm");
-  desc.add<double>("SigmaZ")->setComment("in cm");
-  desc.add<double>("BetaStar")->setComment("in cm");
-  desc.add<double>("Emittance")->setComment("in cm");
+  desc.add<double>("X0", 0.0)->setComment("in cm");
+  desc.add<double>("Y0", 0.0)->setComment("in cm");
+  desc.add<double>("Z0", 0.0)->setComment("in cm");
+  desc.add<double>("MeanX", 0.0)->setComment("in cm");
+  desc.add<double>("MeanY", 0.0)->setComment("in cm");
+  desc.add<double>("MeanZ", 0.0)->setComment("in cm");
+  desc.add<double>("SigmaX", -1.0)->setComment("in cm");
+  desc.add<double>("SigmaY", -1.0)->setComment("in cm");
+  desc.add<double>("SigmaZ", 0.0)->setComment("in cm");
+  desc.add<double>("BetaStar", 0.0)->setComment("in cm");
+  desc.add<double>("Emittance", 0.0)->setComment("in cm");
   desc.add<double>("Alpha", 0.0)->setComment("in radians");
   desc.add<double>("Phi", 0.0)->setComment("in radians");
   desc.add<double>("TimeOffset", 0.0)->setComment("in ns");

CondTools/BeamSpot/test/BeamProfile2DBWriterAll_cfg.py

@@ -31,20 +31,20 @@
 from CondTools.BeamSpot.beamProfile2DBWriter_cfi import beamProfile2DBWriter
 process.load("IOMC.EventVertexGenerators.VtxSmearedParameters_cfi")
 
-# Get the PSets that have BetaStar as a parameter
-psets_with_BetaStar = []
+# Get the PSets that have BetaStar (BetafuncEvtVtxGenerator) or SigmaX (GaussEvtVtxGenerator) as a parameter
+psets_SimBeamSpotObjects = []
 psets_names = []
 for psetName in process.__dict__:
     pset = getattr(process, psetName)
-    if isinstance(pset, cms.PSet) and "BetaStar" in pset.parameterNames_():
+    if isinstance(pset, cms.PSet) and ("BetaStar" in pset.parameterNames_() or "SigmaX" in pset.parameterNames_()):
         print(psetName)
         psets_names.append(psetName)
-        psets_with_BetaStar.append(pset)
+        psets_SimBeamSpotObjects.append(pset)
 
 # Create a VPSet to store the parameter sets
 myVPSet = cms.VPSet()
 
-for i, pset in enumerate(psets_with_BetaStar):
+for i, pset in enumerate(psets_SimBeamSpotObjects):
     cloneName = 'BeamProfile2DBWriter_' + str(i+1)  # Unique clone name
     setattr(process, cloneName, beamProfile2DBWriter.clone(pset,
                                                            recordName = cms.string(psets_names[i])))
@@ -74,7 +74,7 @@
 process.end = cms.EndPath()
 
 process.schedule = cms.Schedule()
-for i, pset in enumerate(psets_with_BetaStar):
+for i, pset in enumerate(psets_SimBeamSpotObjects):
     pathName = 'Path_' + str(i+1)  # Unique path name
     process.schedule.append(getattr(process, pathName))
 process.schedule.append(process.end)

CondTools/BeamSpot/test/BeamProfile2DBWriter_cfg.py

@@ -42,11 +42,22 @@
 process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(1))
 
 from CondTools.BeamSpot.beamProfile2DBWriter_cfi import beamProfile2DBWriter
+# For the BetaFunc smearing (realistic Beamspot) set the following parameters:
+# - X0, Y0, Z0, SigmaZ, BetaStar, Emittance
 process.BeamProfile2DBWriter = beamProfile2DBWriter.clone(X0        = 0.0458532,
                                                           Y0        = -0.016966,
                                                           Z0        = -0.074992,
                                                           SigmaZ    = 3.6,
                                                           BetaStar  = 30.0,
                                                           Emittance = 3.931e-8,)
 
+# For the Gaussian smearing (ideal Beamspot) set the following parameters:
+# - MeanX, MeanY, MeanZ, SigmaX, SigmaY, SigmaZ
+#process.BeamProfile2DBWriter = beamProfile2DBWriter.clone(MeanX     = 0.0,
+#                                                          MeanY     = 0.0,
+#                                                          MeanZ     = 0.0,
+#                                                          SigmaX    = 0.0015,
+#                                                          SigmaY    = 0.0015,
+#                                                          SigmaZ    = 3.6,)
+
 process.p = cms.Path(process.BeamProfile2DBWriter)