From 9cff8571d072757bc554477b3002da405c05f856 Mon Sep 17 00:00:00 2001 From: francescobrivio Date: Fri, 3 Nov 2023 22:50:21 +0100 Subject: [PATCH] update HLLHCEvtVtxGenerator to read BeamSpot from DB --- .../interface/HLLHCEvtVtxGenerator.h | 90 +++++------- .../src/HLLHCEvtVtxGenerator.cc | 134 ++++++++++++------ 2 files changed, 124 insertions(+), 100 deletions(-) diff --git a/IOMC/EventVertexGenerators/interface/HLLHCEvtVtxGenerator.h b/IOMC/EventVertexGenerators/interface/HLLHCEvtVtxGenerator.h index 40f2de9ce5ed0..0ff5befbeb92b 100644 --- a/IOMC/EventVertexGenerators/interface/HLLHCEvtVtxGenerator.h +++ b/IOMC/EventVertexGenerators/interface/HLLHCEvtVtxGenerator.h @@ -12,6 +12,10 @@ */ #include "IOMC/EventVertexGenerators/interface/BaseEvtVtxGenerator.h" +#include "FWCore/Framework/interface/ESWatcher.h" +#include "FWCore/Utilities/interface/ESGetToken.h" +#include "CondFormats/DataRecord/interface/SimBeamSpotHLLHCObjectsRcd.h" +#include "CondFormats/BeamSpotObjects/interface/SimBeamSpotHLLHCObjects.h" #include @@ -33,7 +37,9 @@ class HLLHCEvtVtxGenerator : public BaseEvtVtxGenerator { /** Copy assignment operator */ HLLHCEvtVtxGenerator& operator=(const HLLHCEvtVtxGenerator& rhs) = delete; - ~HLLHCEvtVtxGenerator() override; + ~HLLHCEvtVtxGenerator() override = default; + + void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override; static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); @@ -43,59 +49,29 @@ class HLLHCEvtVtxGenerator : public BaseEvtVtxGenerator { TMatrixD const* GetInvLorentzBoost() const override { return nullptr; }; private: - //spatial and time offset for mean collision - const double fMeanX, fMeanY, fMeanZ, fTimeOffset; - - //proton beam energy - const double momeV; - const double gamma; - const double beta; - const double betagamma; - - //crossing angle - const double phi; - - //crab cavity frequency - const double wcc; - - // 800 MHz RF? - const bool RF800; - - //beta crossing plane (m) - const double betx; - - //beta separation plane (m) - const double bets; - - //horizontal emittance - const double epsxn; - - //vertical emittance - const double epssn; - - //bunch length - const double sigs; - - //crabbing angle crossing - const double alphax; - - //crabbing angle separation - const double alphay; - - // ratio of crabbing angle to crossing angle - const double oncc; - - //normalized crossing emittance - const double epsx; - - //normlaized separation emittance - const double epss; - - //size in x - const double sigx; - - // crossing angle * crab frequency - const double phiCR; + // Configurable parameters + double fMeanX, fMeanY, fMeanZ, fTimeOffset; //spatial and time offset for mean collision + double fEProton; // proton beam energy + double fCrossingAngle; // crossing angle + double fCrabFrequency; // crab cavity frequency + bool fRF800; // 800 MHz RF? + double fBetaCrossingPlane; // beta crossing plane (m) + double fBetaSeparationPlane; // beta separation plane (m) + double fHorizontalEmittance; // horizontal emittance + double fVerticalEmittance; // vertical emittance + double fBunchLength; // bunch length + double fCrabbingAngleCrossing; // crabbing angle crossing + double fCrabbingAngleSeparation; // crabbing angle separation + + // Parameters inferred from configurables + double gamma; // beam configurations + double beta; + double betagamma; + double oncc; // ratio of crabbing angle to crossing angle + double epsx; // normalized crossing emittance + double epss; // normlaized separation emittance + double sigx; // size in x + double phiCR; // crossing angle * crab frequency //width for y plane double sigma(double z, double epsilon, double beta, double betagamma) const; @@ -105,6 +81,12 @@ class HLLHCEvtVtxGenerator : public BaseEvtVtxGenerator { // 4D intensity double intensity(double x, double y, double z, double t) const; + + // Read frm DB + bool readDB_; + void update(const edm::EventSetup& iEventSetup); + edm::ESWatcher parameterWatcher_; + edm::ESGetToken beamToken_; }; #endif diff --git a/IOMC/EventVertexGenerators/src/HLLHCEvtVtxGenerator.cc b/IOMC/EventVertexGenerators/src/HLLHCEvtVtxGenerator.cc index 48ffec65f909a..ea4a56e4e064b 100644 --- a/IOMC/EventVertexGenerators/src/HLLHCEvtVtxGenerator.cc +++ b/IOMC/EventVertexGenerators/src/HLLHCEvtVtxGenerator.cc @@ -16,7 +16,6 @@ using namespace std; namespace { - constexpr double pmass = 0.9382720813e9; // eV constexpr double gamma34 = 1.22541670246517764513; // Gamma(3/4) constexpr double gamma14 = 3.62560990822190831193; // Gamma(1/4) @@ -48,35 +47,77 @@ void HLLHCEvtVtxGenerator::fillDescriptions(edm::ConfigurationDescriptions& desc descriptions.add("HLLHCEvtVtxGenerator", desc); } -HLLHCEvtVtxGenerator::HLLHCEvtVtxGenerator(const edm::ParameterSet& p) - : BaseEvtVtxGenerator(p), - fMeanX(p.getParameter("MeanXIncm") * cm), - fMeanY(p.getParameter("MeanYIncm") * cm), - fMeanZ(p.getParameter("MeanZIncm") * cm), - fTimeOffset(p.getParameter("TimeOffsetInns") * ns * c_light), - momeV(p.getParameter("EprotonInGeV") * 1e9), - gamma(momeV / pmass + 1.0), - beta(std::sqrt((1.0 - 1.0 / gamma) * ((1.0 + 1.0 / gamma)))), - betagamma(beta * gamma), - phi(p.getParameter("CrossingAngleInurad") * 1e-6), - wcc(p.getParameter("CrabFrequencyInMHz") * 1e6), - RF800(p.getParameter("RF800")), - betx(p.getParameter("BetaCrossingPlaneInm")), - bets(p.getParameter("BetaSeparationPlaneInm")), - epsxn(p.getParameter("HorizontalEmittance")), - epssn(p.getParameter("VerticalEmittance")), - sigs(p.getParameter("BunchLengthInm")), - alphax(p.getParameter("CrabbingAngleCrossingInurad") * 1e-6), - alphay(p.getParameter("CrabbingAngleSeparationInurad") * 1e-6), - oncc(alphax / phi), - epsx(epsxn / (betagamma)), - epss(epsx), - sigx(std::sqrt(epsx * betx)), - phiCR(oncc * phi) - -{} - -HLLHCEvtVtxGenerator::~HLLHCEvtVtxGenerator() {} +HLLHCEvtVtxGenerator::HLLHCEvtVtxGenerator(const edm::ParameterSet& p) : BaseEvtVtxGenerator(p) { + readDB_ = p.getParameter("readDB"); + if (!readDB_) { + // Read configurable parameters + fMeanX = p.getParameter("MeanXIncm") * cm; + fMeanY = p.getParameter("MeanYIncm") * cm; + fMeanZ = p.getParameter("MeanZIncm") * cm; + fTimeOffset = p.getParameter("TimeOffsetInns") * ns * c_light; + fEProton = p.getParameter("EprotonInGeV") * 1e9; + fCrossingAngle = p.getParameter("CrossingAngleInurad") * 1e-6; + fCrabFrequency = p.getParameter("CrabFrequencyInMHz") * 1e6; + fRF800 = p.getParameter("RF800"); + fBetaCrossingPlane = p.getParameter("BetaCrossingPlaneInm"); + fBetaSeparationPlane = p.getParameter("BetaSeparationPlaneInm"); + fHorizontalEmittance = p.getParameter("HorizontalEmittance"); + fVerticalEmittance = p.getParameter("VerticalEmittance"); + fBunchLength = p.getParameter("BunchLengthInm"); + fCrabbingAngleCrossing = p.getParameter("CrabbingAngleCrossingInurad") * 1e-6; + fCrabbingAngleSeparation = p.getParameter("CrabbingAngleSeparationInurad") * 1e-6; + // Set parameters inferred from configurables + gamma = fEProton / pmass + 1.0; + beta = std::sqrt((1.0 - 1.0 / gamma) * ((1.0 + 1.0 / gamma))); + betagamma = beta * gamma; + oncc = fCrabbingAngleCrossing / fCrossingAngle; + epsx = fHorizontalEmittance / (betagamma); + epss = epsx; + sigx = std::sqrt(epsx * fBetaCrossingPlane); + phiCR = oncc * fCrossingAngle; + } + if (readDB_) { + // NOTE: this is currently watching LS transitions, while it should watch Run transitions, + // even though in reality there is no Run Dependent MC (yet) in CMS + beamToken_ = + esConsumes(); + } +} + +void HLLHCEvtVtxGenerator::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const& iEventSetup) { + update(iEventSetup); +} + +void HLLHCEvtVtxGenerator::update(const edm::EventSetup& iEventSetup) { + if (readDB_ && parameterWatcher_.check(iEventSetup)) { + edm::ESHandle beamhandle = iEventSetup.getHandle(beamToken_); + // Read configurable parameters + fMeanX = beamhandle->meanX() * cm; + fMeanY = beamhandle->meanY() * cm; + fMeanZ = beamhandle->meanZ() * cm; + fEProton = beamhandle->eProton() * 1e9; + fCrossingAngle = beamhandle->crossingAngle() * 1e-6; + fCrabFrequency = beamhandle->crabFrequency() * 1e6; + fRF800 = beamhandle->rf800(); + fBetaCrossingPlane = beamhandle->betaCrossingPlane(); + fBetaSeparationPlane = beamhandle->betaSeparationPlane(); + fHorizontalEmittance = beamhandle->horizontalEmittance(); + fVerticalEmittance = beamhandle->verticalEmittance(); + fBunchLength = beamhandle->bunchLenght(); + fCrabbingAngleCrossing = beamhandle->crabbingAngleCrossing() * 1e-6; + fCrabbingAngleSeparation = beamhandle->crabbingAngleSeparation() * 1e-6; + fTimeOffset = beamhandle->timeOffset() * ns * c_light; + // Set parameters inferred from configurables + gamma = fEProton / pmass + 1.0; + beta = std::sqrt((1.0 - 1.0 / gamma) * ((1.0 + 1.0 / gamma))); + betagamma = beta * gamma; + oncc = fCrabbingAngleCrossing / fCrossingAngle; + epsx = fHorizontalEmittance / (betagamma); + epss = epsx; + sigx = std::sqrt(epsx * fBetaCrossingPlane); + phiCR = oncc * fCrossingAngle; + } +} HepMC::FourVector HLLHCEvtVtxGenerator::newVertex(CLHEP::HepRandomEngine* engine) const { double imax = intensity(0., 0., 0., 0.); @@ -88,10 +129,10 @@ HepMC::FourVector HLLHCEvtVtxGenerator::newVertex(CLHEP::HepRandomEngine* engine auto shoot = [&]() { return CLHEP::RandFlat::shoot(engine); }; do { - z = (shoot() - 0.5) * 6.0 * sigs; - t = (shoot() - 0.5) * 6.0 * sigs; - x = (shoot() - 0.5) * 12.0 * sigma(0.0, epsxn, betx, betagamma); - y = (shoot() - 0.5) * 12.0 * sigma(0.0, epssn, bets, betagamma); + z = (shoot() - 0.5) * 6.0 * fBunchLength; + t = (shoot() - 0.5) * 6.0 * fBunchLength; + x = (shoot() - 0.5) * 12.0 * sigma(0.0, fHorizontalEmittance, fBetaCrossingPlane, betagamma); + y = (shoot() - 0.5) * 12.0 * sigma(0.0, fVerticalEmittance, fBetaSeparationPlane, betagamma); i = intensity(x, y, z, t); @@ -119,7 +160,6 @@ HepMC::FourVector HLLHCEvtVtxGenerator::newVertex(CLHEP::HepRandomEngine* engine double HLLHCEvtVtxGenerator::sigma(double z, double epsilon, double beta, double betagamma) const { double sigma = std::sqrt(epsilon * (beta + z * z / beta) / betagamma); - return sigma; } @@ -127,9 +167,9 @@ double HLLHCEvtVtxGenerator::intensity(double x, double y, double z, double t) c //---c in m/s --- remember t is already in meters constexpr double c = 2.99792458e+8; // m/s - const double sigmay = sigma(z, epssn, bets, betagamma); + const double sigmay = sigma(z, fVerticalEmittance, fBetaSeparationPlane, betagamma); - const double alphay_mod = alphay * std::cos(wcc * (z - t) / c); + const double alphay_mod = fCrabbingAngleSeparation * std::cos(fCrabFrequency * (z - t) / c); const double cay = std::cos(alphay_mod); const double say = std::sin(alphay_mod); @@ -146,17 +186,17 @@ double HLLHCEvtVtxGenerator::intensity(double x, double y, double z, double t) c double HLLHCEvtVtxGenerator::integrandCC(double x, double z, double ct) const { constexpr double local_c_light = 2.99792458e8; - const double k = wcc / local_c_light * two_pi; + const double k = fCrabFrequency / local_c_light * two_pi; const double k2 = k * k; - const double cos = std::cos(phi / 2.0); - const double sin = std::sin(phi / 2.0); + const double cos = std::cos(fCrossingAngle / 2.0); + const double sin = std::sin(fCrossingAngle / 2.0); const double cos2 = cos * cos; const double sin2 = sin * sin; const double sigx2 = sigx * sigx; - const double sigmax2 = sigx2 * (1 + z * z / (betx * betx)); + const double sigmax2 = sigx2 * (1 + z * z / (fBetaCrossingPlane * fBetaCrossingPlane)); - const double sigs2 = sigs * sigs; + const double sigs2 = fBunchLength * fBunchLength; constexpr double factorRMSgauss4 = 1. / sqrt2 / gamma34 * gamma14; // # Factor to take rms sigma as input of the supergaussian @@ -167,7 +207,7 @@ double HLLHCEvtVtxGenerator::integrandCC(double x, double z, double ct) const { double result = -1.0; - if (!RF800) { + if (!fRF800) { const double norm = 2.0 / (two_pi * sigs2); const double cosks = std::cos(k * z); const double sinkct = std::sin(k * ct); @@ -183,9 +223,10 @@ double HLLHCEvtVtxGenerator::integrandCC(double x, double z, double ct) const { + x * ct * sin / sigs2 // contribution from x integrand + 2 * x * cos * cosks * sinkct * sinCR / k / sigmax2 // contribution from x integrand //+(2*ct/k)*np.cos(k*s)*np.sin(k*ct) *(sin*sinCR)/(sigs2*cos) # small term - //+ct**2*(sin2/sigs4)/(cos2/sigmax2) # small term + //+ct**2*(sin2/sigs4)/(cos2/sigmax2) # small term ) / - (1.0 + (z * z) / (betx * betx)) / std::sqrt(1.0 + (z * z) / (bets * bets)); + (1.0 + (z * z) / (fBetaCrossingPlane * fBetaCrossingPlane)) / + std::sqrt(1.0 + (z * z) / (fBetaSeparationPlane * fBetaSeparationPlane)); } else { const double norm = 2.0 / (NormFactorGauss4 * sigs2 * factorRMSgauss4); @@ -198,7 +239,8 @@ double HLLHCEvtVtxGenerator::integrandCC(double x, double z, double ct) const { sin2 / (4 * k2 * sigmax2) * (2 + 4 * k2 * z * z - std::cos(2 * k * (z - ct)) - std::cos(2 * k * (z + ct)) - 8 * k * s * std::cos(k * ct) * std::sin(k * z) - 4 * cosks * cosks * sinct * sinct)) / - std::sqrt(1 + z * z / (betx * betx)) / std::sqrt(1 + z * z / (bets * bets)); + std::sqrt(1 + z * z / (fBetaCrossingPlane * fBetaCrossingPlane)) / + std::sqrt(1 + z * z / (fBetaSeparationPlane * fBetaSeparationPlane)); } return result;