Skip to content

Commit

Permalink
update HLLHCEvtVtxGenerator to read BeamSpot from DB
Browse files Browse the repository at this point in the history
  • Loading branch information
@francescobrivio
francescobrivio committed Nov 3, 2023
1 parent 26730cd commit 9cff857
Show file tree
Hide file tree
Showing 2 changed files with 124 additions and 100 deletions.
90 changes: 36 additions & 54 deletions IOMC/EventVertexGenerators/interface/HLLHCEvtVtxGenerator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <string>

Expand All @@ -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);

Expand All @@ -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;
Expand All @@ -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<SimBeamSpotHLLHCObjectsRcd> parameterWatcher_;
edm::ESGetToken<SimBeamSpotHLLHCObjects, SimBeamSpotHLLHCObjectsRcd> beamToken_;
};

#endif
134 changes: 88 additions & 46 deletions IOMC/EventVertexGenerators/src/HLLHCEvtVtxGenerator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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<double>("MeanXIncm") * cm),
fMeanY(p.getParameter<double>("MeanYIncm") * cm),
fMeanZ(p.getParameter<double>("MeanZIncm") * cm),
fTimeOffset(p.getParameter<double>("TimeOffsetInns") * ns * c_light),
momeV(p.getParameter<double>("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<double>("CrossingAngleInurad") * 1e-6),
wcc(p.getParameter<double>("CrabFrequencyInMHz") * 1e6),
RF800(p.getParameter<bool>("RF800")),
betx(p.getParameter<double>("BetaCrossingPlaneInm")),
bets(p.getParameter<double>("BetaSeparationPlaneInm")),
epsxn(p.getParameter<double>("HorizontalEmittance")),
epssn(p.getParameter<double>("VerticalEmittance")),
sigs(p.getParameter<double>("BunchLengthInm")),
alphax(p.getParameter<double>("CrabbingAngleCrossingInurad") * 1e-6),
alphay(p.getParameter<double>("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<bool>("readDB");
if (!readDB_) {
// Read configurable parameters
fMeanX = p.getParameter<double>("MeanXIncm") * cm;
fMeanY = p.getParameter<double>("MeanYIncm") * cm;
fMeanZ = p.getParameter<double>("MeanZIncm") * cm;
fTimeOffset = p.getParameter<double>("TimeOffsetInns") * ns * c_light;
fEProton = p.getParameter<double>("EprotonInGeV") * 1e9;
fCrossingAngle = p.getParameter<double>("CrossingAngleInurad") * 1e-6;
fCrabFrequency = p.getParameter<double>("CrabFrequencyInMHz") * 1e6;
fRF800 = p.getParameter<bool>("RF800");
fBetaCrossingPlane = p.getParameter<double>("BetaCrossingPlaneInm");
fBetaSeparationPlane = p.getParameter<double>("BetaSeparationPlaneInm");
fHorizontalEmittance = p.getParameter<double>("HorizontalEmittance");
fVerticalEmittance = p.getParameter<double>("VerticalEmittance");
fBunchLength = p.getParameter<double>("BunchLengthInm");
fCrabbingAngleCrossing = p.getParameter<double>("CrabbingAngleCrossingInurad") * 1e-6;
fCrabbingAngleSeparation = p.getParameter<double>("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<SimBeamSpotHLLHCObjects, SimBeamSpotHLLHCObjectsRcd, edm::Transition::BeginLuminosityBlock>();
}
}

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<SimBeamSpotHLLHCObjects> 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.);
Expand All @@ -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);

Expand Down Expand Up @@ -119,17 +160,16 @@ 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;
}

double HLLHCEvtVtxGenerator::intensity(double x, double y, double z, double t) const {
//---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);
Expand All @@ -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
Expand All @@ -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);
Expand All @@ -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);
Expand All @@ -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;
Expand Down

0 comments on commit 9cff857

Please sign in to comment.