File indexing completed on 2024-12-10 02:31:54
0001
0002
0003 #include "IOMC/EventVertexGenerators/interface/HLLHCEvtVtxGenerator.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include <CLHEP/Random/RandFlat.h>
0012 #include <CLHEP/Units/SystemOfUnits.h>
0013 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0014
0015 using namespace std;
0016
0017 namespace {
0018 constexpr double pmass = 0.9382720813e9;
0019 constexpr double gamma34 = 1.22541670246517764513;
0020 constexpr double gamma14 = 3.62560990822190831193;
0021 constexpr double gamma54 = 0.90640247705547798267;
0022 constexpr double sqrt2 = 1.41421356237;
0023 constexpr double sqrt2to5 = 5.65685424949;
0024 constexpr double two_pi = 2.0 * M_PI;
0025 }
0026
0027 void HLLHCEvtVtxGenerator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0028 edm::ParameterSetDescription desc;
0029 desc.add<double>("MeanXIncm", 0.0);
0030 desc.add<double>("MeanYIncm", 0.0);
0031 desc.add<double>("MeanZIncm", 0.0);
0032 desc.add<double>("TimeOffsetInns", 0.0);
0033 desc.add<double>("EprotonInGeV", 7000.0);
0034 desc.add<double>("CrossingAngleInurad", 510.0);
0035 desc.add<double>("CrabbingAngleCrossingInurad", 380.0);
0036 desc.add<double>("CrabbingAngleSeparationInurad", 0.0);
0037 desc.add<double>("CrabFrequencyInMHz", 400.0);
0038 desc.add<bool>("RF800", false);
0039 desc.add<double>("BetaCrossingPlaneInm", 0.20);
0040 desc.add<double>("BetaSeparationPlaneInm", 0.20);
0041 desc.add<double>("HorizontalEmittance", 2.5e-06);
0042 desc.add<double>("VerticalEmittance", 2.05e-06);
0043 desc.add<double>("BunchLengthInm", 0.09);
0044 desc.add<edm::InputTag>("src");
0045 desc.add<bool>("readDB");
0046 descriptions.add("HLLHCEvtVtxGenerator", desc);
0047 }
0048
0049 HLLHCEvtVtxGenerator::HLLHCEvtVtxGenerator(const edm::ParameterSet& p) : BaseEvtVtxGenerator(p) {
0050 readDB_ = p.getParameter<bool>("readDB");
0051 if (!readDB_) {
0052
0053 fMeanX = p.getParameter<double>("MeanXIncm") * CLHEP::cm;
0054 fMeanY = p.getParameter<double>("MeanYIncm") * CLHEP::cm;
0055 fMeanZ = p.getParameter<double>("MeanZIncm") * CLHEP::cm;
0056 fTimeOffset_c_light = p.getParameter<double>("TimeOffsetInns") * CLHEP::ns * CLHEP::c_light;
0057 fEProton = p.getParameter<double>("EprotonInGeV") * 1e9;
0058 fCrossingAngle = p.getParameter<double>("CrossingAngleInurad") * 1e-6;
0059 fCrabFrequency = p.getParameter<double>("CrabFrequencyInMHz") * 1e6;
0060 fRF800 = p.getParameter<bool>("RF800");
0061 fBetaCrossingPlane = p.getParameter<double>("BetaCrossingPlaneInm");
0062 fBetaSeparationPlane = p.getParameter<double>("BetaSeparationPlaneInm");
0063 fHorizontalEmittance = p.getParameter<double>("HorizontalEmittance");
0064 fVerticalEmittance = p.getParameter<double>("VerticalEmittance");
0065 fBunchLength = p.getParameter<double>("BunchLengthInm");
0066 fCrabbingAngleCrossing = p.getParameter<double>("CrabbingAngleCrossingInurad") * 1e-6;
0067 fCrabbingAngleSeparation = p.getParameter<double>("CrabbingAngleSeparationInurad") * 1e-6;
0068
0069 gamma = fEProton / pmass;
0070 beta = std::sqrt((1.0 - 1.0 / gamma) * ((1.0 + 1.0 / gamma)));
0071 betagamma = beta * gamma;
0072 oncc = fCrabbingAngleCrossing / fCrossingAngle;
0073 epsx = fHorizontalEmittance / (betagamma);
0074 epss = epsx;
0075 sigx = std::sqrt(epsx * fBetaCrossingPlane);
0076 phiCR = oncc * fCrossingAngle;
0077 }
0078 if (readDB_) {
0079
0080
0081 beamToken_ =
0082 esConsumes<SimBeamSpotHLLHCObjects, SimBeamSpotHLLHCObjectsRcd, edm::Transition::BeginLuminosityBlock>();
0083 }
0084 }
0085
0086 void HLLHCEvtVtxGenerator::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const& iEventSetup) {
0087 update(iEventSetup);
0088 }
0089
0090 void HLLHCEvtVtxGenerator::update(const edm::EventSetup& iEventSetup) {
0091 if (readDB_ && parameterWatcher_.check(iEventSetup)) {
0092 edm::ESHandle<SimBeamSpotHLLHCObjects> beamhandle = iEventSetup.getHandle(beamToken_);
0093
0094 fMeanX = beamhandle->meanX() * CLHEP::cm;
0095 fMeanY = beamhandle->meanY() * CLHEP::cm;
0096 fMeanZ = beamhandle->meanZ() * CLHEP::cm;
0097 fEProton = beamhandle->eProton() * 1e9;
0098 fCrossingAngle = beamhandle->crossingAngle() * 1e-6;
0099 fCrabFrequency = beamhandle->crabFrequency() * 1e6;
0100 fRF800 = beamhandle->rf800();
0101 fBetaCrossingPlane = beamhandle->betaCrossingPlane();
0102 fBetaSeparationPlane = beamhandle->betaSeparationPlane();
0103 fHorizontalEmittance = beamhandle->horizontalEmittance();
0104 fVerticalEmittance = beamhandle->verticalEmittance();
0105 fBunchLength = beamhandle->bunchLenght();
0106 fCrabbingAngleCrossing = beamhandle->crabbingAngleCrossing() * 1e-6;
0107 fCrabbingAngleSeparation = beamhandle->crabbingAngleSeparation() * 1e-6;
0108 fTimeOffset_c_light = beamhandle->timeOffset() * CLHEP::ns * CLHEP::c_light;
0109
0110 gamma = fEProton / pmass;
0111 beta = std::sqrt((1.0 - 1.0 / gamma) * ((1.0 + 1.0 / gamma)));
0112 betagamma = beta * gamma;
0113 oncc = fCrabbingAngleCrossing / fCrossingAngle;
0114 epsx = fHorizontalEmittance / (betagamma);
0115 epss = epsx;
0116 sigx = std::sqrt(epsx * fBetaCrossingPlane);
0117 phiCR = oncc * fCrossingAngle;
0118 }
0119 }
0120
0121 ROOT::Math::XYZTVector HLLHCEvtVtxGenerator::vertexShift(CLHEP::HepRandomEngine* engine) const {
0122 double imax = intensity(0., 0., 0., 0.);
0123
0124 double x(0.), y(0.), z(0.), t(0.), i(0.);
0125
0126 int count = 0;
0127
0128 auto shoot = [&]() { return CLHEP::RandFlat::shoot(engine); };
0129
0130 do {
0131 z = (shoot() - 0.5) * 6.0 * fBunchLength;
0132 t = (shoot() - 0.5) * 6.0 * fBunchLength;
0133 x = (shoot() - 0.5) * 12.0 * sigma(0.0, fHorizontalEmittance, fBetaCrossingPlane, betagamma);
0134 y = (shoot() - 0.5) * 12.0 * sigma(0.0, fVerticalEmittance, fBetaSeparationPlane, betagamma);
0135
0136 i = intensity(x, y, z, t);
0137
0138 if (i > imax)
0139 edm::LogError("Too large intensity") << "i>imax : " << i << " > " << imax << endl;
0140 ++count;
0141 } while ((i < imax * shoot()) && count < 10000);
0142
0143 if (count > 9999)
0144 edm::LogError("Too many tries ") << " count : " << count << endl;
0145
0146
0147 x *= 1000.0;
0148 y *= 1000.0;
0149 z *= 1000.0;
0150 t *= 1000.0;
0151
0152 x += fMeanX;
0153 y += fMeanY;
0154 z += fMeanZ;
0155 t += fTimeOffset_c_light;
0156
0157 return ROOT::Math::XYZTVector(x, y, z, t);
0158 }
0159
0160 double HLLHCEvtVtxGenerator::sigma(double z, double epsilon, double beta, double betagamma) const {
0161 double sigma = std::sqrt(epsilon * (beta + z * z / beta) / betagamma);
0162 return sigma;
0163 }
0164
0165 double HLLHCEvtVtxGenerator::intensity(double x, double y, double z, double t) const {
0166
0167 constexpr double c = 2.99792458e+8;
0168
0169 const double sigmay = sigma(z, fVerticalEmittance, fBetaSeparationPlane, betagamma);
0170
0171 const double alphay_mod = fCrabbingAngleSeparation * std::cos(fCrabFrequency * (z - t) / c);
0172
0173 const double cay = std::cos(alphay_mod);
0174 const double say = std::sin(alphay_mod);
0175
0176 const double dy = -(z - t) * say - y * cay;
0177
0178 const double xzt_density = integrandCC(x, z, t);
0179
0180 const double norm = two_pi * sigmay;
0181
0182 return (std::exp(-dy * dy / (sigmay * sigmay)) * xzt_density / norm);
0183 }
0184
0185 double HLLHCEvtVtxGenerator::integrandCC(double x, double z, double ct) const {
0186 constexpr double local_c_light = 2.99792458e8;
0187
0188 const double k = fCrabFrequency / local_c_light * two_pi;
0189 const double k2 = k * k;
0190 const double cos = std::cos(fCrossingAngle / 2.0);
0191 const double sin = std::sin(fCrossingAngle / 2.0);
0192 const double cos2 = cos * cos;
0193 const double sin2 = sin * sin;
0194
0195 const double sigx2 = sigx * sigx;
0196 const double sigmax2 = sigx2 * (1 + z * z / (fBetaCrossingPlane * fBetaCrossingPlane));
0197
0198 const double sigs2 = fBunchLength * fBunchLength;
0199
0200 constexpr double factorRMSgauss4 =
0201 1. / sqrt2 / gamma34 * gamma14;
0202 constexpr double NormFactorGauss4 = sqrt2to5 * gamma54 * gamma54;
0203
0204 const double sinCR = std::sin(phiCR / 2.0);
0205 const double sinCR2 = sinCR * sinCR;
0206
0207 double result = -1.0;
0208
0209 if (!fRF800) {
0210 const double norm = 2.0 / (two_pi * sigs2);
0211 const double cosks = std::cos(k * z);
0212 const double sinkct = std::sin(k * ct);
0213 result = norm *
0214 std::exp(-ct * ct / sigs2 - z * z * cos2 / sigs2 -
0215 1.0 / (4 * k2 * sigmax2) *
0216 (
0217
0218 -8 * z * k * std::sin(k * z) * std::cos(k * ct) * sin * sinCR + 2 * sinCR2 -
0219 std::cos(2 * k * (z - ct)) * sinCR2 - std::cos(2 * k * (z + ct)) * sinCR2 +
0220 4 * k2 * z * z * sin2) -
0221 x * x * (cos2 / sigmax2 + sin2 / sigs2)
0222 + x * ct * sin / sigs2
0223 + 2 * x * cos * cosks * sinkct * sinCR / k / sigmax2
0224
0225
0226 ) /
0227 (1.0 + (z * z) / (fBetaCrossingPlane * fBetaCrossingPlane)) /
0228 std::sqrt(1.0 + (z * z) / (fBetaSeparationPlane * fBetaSeparationPlane));
0229
0230 } else {
0231 const double norm = 2.0 / (NormFactorGauss4 * sigs2 * factorRMSgauss4);
0232 const double sigs4 = sigs2 * sigs2 * factorRMSgauss4 * factorRMSgauss4;
0233 const double cosks = std::cos(k * z);
0234 const double sinct = std::sin(k * ct);
0235 result =
0236 norm *
0237 std::exp(-ct * ct * ct * ct / sigs4 - z * z * z * z * cos2 * cos2 / sigs4 - 6 * ct * ct * z * z * cos2 / sigs4 -
0238 sin2 / (4 * k2 * sigmax2) *
0239 (2 + 4 * k2 * z * z - std::cos(2 * k * (z - ct)) - std::cos(2 * k * (z + ct)) -
0240 8 * k * CLHEP::s * std::cos(k * ct) * std::sin(k * z) - 4 * cosks * cosks * sinct * sinct)) /
0241 std::sqrt((1 + z * z / (fBetaCrossingPlane * fBetaCrossingPlane)) /
0242 (1 + z * z / (fBetaSeparationPlane * fBetaSeparationPlane)));
0243 }
0244
0245 return result;
0246 }