File indexing completed on 2023-03-17 11:10:13
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/GlobalSystemOfUnits.h"
0013 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0014 #include "HepMC/SimpleVector.h"
0015
0016 using namespace std;
0017
0018 namespace {
0019
0020 constexpr double pmass = 0.9382720813e9;
0021 constexpr double gamma34 = 1.22541670246517764513;
0022 constexpr double gamma14 = 3.62560990822190831193;
0023 constexpr double gamma54 = 0.90640247705547798267;
0024 constexpr double sqrt2 = 1.41421356237;
0025 constexpr double sqrt2to5 = 5.65685424949;
0026 constexpr double two_pi = 2.0 * M_PI;
0027 }
0028
0029 void HLLHCEvtVtxGenerator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0030 edm::ParameterSetDescription desc;
0031 desc.add<double>("MeanXIncm", 0.0);
0032 desc.add<double>("MeanYIncm", 0.0);
0033 desc.add<double>("MeanZIncm", 0.0);
0034 desc.add<double>("TimeOffsetInns", 0.0);
0035 desc.add<double>("EprotonInGeV", 7000.0);
0036 desc.add<double>("CrossingAngleInurad", 510.0);
0037 desc.add<double>("CrabbingAngleCrossingInurad", 380.0);
0038 desc.add<double>("CrabbingAngleSeparationInurad", 0.0);
0039 desc.add<double>("CrabFrequencyInMHz", 400.0);
0040 desc.add<bool>("RF800", false);
0041 desc.add<double>("BetaCrossingPlaneInm", 0.20);
0042 desc.add<double>("BetaSeparationPlaneInm", 0.20);
0043 desc.add<double>("HorizontalEmittance", 2.5e-06);
0044 desc.add<double>("VerticalEmittance", 2.05e-06);
0045 desc.add<double>("BunchLengthInm", 0.09);
0046 desc.add<edm::InputTag>("src");
0047 desc.add<bool>("readDB");
0048 descriptions.add("HLLHCEvtVtxGenerator", desc);
0049 }
0050
0051 HLLHCEvtVtxGenerator::HLLHCEvtVtxGenerator(const edm::ParameterSet& p)
0052 : BaseEvtVtxGenerator(p),
0053 fMeanX(p.getParameter<double>("MeanXIncm") * cm),
0054 fMeanY(p.getParameter<double>("MeanYIncm") * cm),
0055 fMeanZ(p.getParameter<double>("MeanZIncm") * cm),
0056 fTimeOffset(p.getParameter<double>("TimeOffsetInns") * ns * c_light),
0057 momeV(p.getParameter<double>("EprotonInGeV") * 1e9),
0058 gamma(momeV / pmass + 1.0),
0059 beta(std::sqrt((1.0 - 1.0 / gamma) * ((1.0 + 1.0 / gamma)))),
0060 betagamma(beta * gamma),
0061 phi(p.getParameter<double>("CrossingAngleInurad") * 1e-6),
0062 wcc(p.getParameter<double>("CrabFrequencyInMHz") * 1e6),
0063 RF800(p.getParameter<bool>("RF800")),
0064 betx(p.getParameter<double>("BetaCrossingPlaneInm")),
0065 bets(p.getParameter<double>("BetaSeparationPlaneInm")),
0066 epsxn(p.getParameter<double>("HorizontalEmittance")),
0067 epssn(p.getParameter<double>("VerticalEmittance")),
0068 sigs(p.getParameter<double>("BunchLengthInm")),
0069 alphax(p.getParameter<double>("CrabbingAngleCrossingInurad") * 1e-6),
0070 alphay(p.getParameter<double>("CrabbingAngleSeparationInurad") * 1e-6),
0071 oncc(alphax / phi),
0072 epsx(epsxn / (betagamma)),
0073 epss(epsx),
0074 sigx(std::sqrt(epsx * betx)),
0075 phiCR(oncc * phi)
0076
0077 {}
0078
0079 HLLHCEvtVtxGenerator::~HLLHCEvtVtxGenerator() {}
0080
0081 HepMC::FourVector HLLHCEvtVtxGenerator::newVertex(CLHEP::HepRandomEngine* engine) const {
0082 double imax = intensity(0., 0., 0., 0.);
0083
0084 double x(0.), y(0.), z(0.), t(0.), i(0.);
0085
0086 int count = 0;
0087
0088 auto shoot = [&]() { return CLHEP::RandFlat::shoot(engine); };
0089
0090 do {
0091 z = (shoot() - 0.5) * 6.0 * sigs;
0092 t = (shoot() - 0.5) * 6.0 * sigs;
0093 x = (shoot() - 0.5) * 12.0 * sigma(0.0, epsxn, betx, betagamma);
0094 y = (shoot() - 0.5) * 12.0 * sigma(0.0, epssn, bets, betagamma);
0095
0096 i = intensity(x, y, z, t);
0097
0098 if (i > imax)
0099 edm::LogError("Too large intensity") << "i>imax : " << i << " > " << imax << endl;
0100 ++count;
0101 } while ((i < imax * shoot()) && count < 10000);
0102
0103 if (count > 9999)
0104 edm::LogError("Too many tries ") << " count : " << count << endl;
0105
0106
0107 x *= 1000.0;
0108 y *= 1000.0;
0109 z *= 1000.0;
0110 t *= 1000.0;
0111
0112 x += fMeanX;
0113 y += fMeanY;
0114 z += fMeanZ;
0115 t += fTimeOffset;
0116
0117 return HepMC::FourVector(x, y, z, t);
0118 }
0119
0120 double HLLHCEvtVtxGenerator::sigma(double z, double epsilon, double beta, double betagamma) const {
0121 double sigma = std::sqrt(epsilon * (beta + z * z / beta) / betagamma);
0122
0123 return sigma;
0124 }
0125
0126 double HLLHCEvtVtxGenerator::intensity(double x, double y, double z, double t) const {
0127
0128 constexpr double c = 2.99792458e+8;
0129
0130 const double sigmay = sigma(z, epssn, bets, betagamma);
0131
0132 const double alphay_mod = alphay * std::cos(wcc * (z - t) / c);
0133
0134 const double cay = std::cos(alphay_mod);
0135 const double say = std::sin(alphay_mod);
0136
0137 const double dy = -(z - t) * say - y * cay;
0138
0139 const double xzt_density = integrandCC(x, z, t);
0140
0141 const double norm = two_pi * sigmay;
0142
0143 return (std::exp(-dy * dy / (sigmay * sigmay)) * xzt_density / norm);
0144 }
0145
0146 double HLLHCEvtVtxGenerator::integrandCC(double x, double z, double ct) const {
0147 constexpr double local_c_light = 2.99792458e8;
0148
0149 const double k = wcc / local_c_light * two_pi;
0150 const double k2 = k * k;
0151 const double cos = std::cos(phi / 2.0);
0152 const double sin = std::sin(phi / 2.0);
0153 const double cos2 = cos * cos;
0154 const double sin2 = sin * sin;
0155
0156 const double sigx2 = sigx * sigx;
0157 const double sigmax2 = sigx2 * (1 + z * z / (betx * betx));
0158
0159 const double sigs2 = sigs * sigs;
0160
0161 constexpr double factorRMSgauss4 =
0162 1. / sqrt2 / gamma34 * gamma14;
0163 constexpr double NormFactorGauss4 = sqrt2to5 * gamma54 * gamma54;
0164
0165 const double sinCR = std::sin(phiCR / 2.0);
0166 const double sinCR2 = sinCR * sinCR;
0167
0168 double result = -1.0;
0169
0170 if (!RF800) {
0171 const double norm = 2.0 / (two_pi * sigs2);
0172 const double cosks = std::cos(k * z);
0173 const double sinkct = std::sin(k * ct);
0174 result = norm *
0175 std::exp(-ct * ct / sigs2 - z * z * cos2 / sigs2 -
0176 1.0 / (4 * k2 * sigmax2) *
0177 (
0178
0179 -8 * z * k * std::sin(k * z) * std::cos(k * ct) * sin * sinCR + 2 * sinCR2 -
0180 std::cos(2 * k * (z - ct)) * sinCR2 - std::cos(2 * k * (z + ct)) * sinCR2 +
0181 4 * k2 * z * z * sin2) -
0182 x * x * (cos2 / sigmax2 + sin2 / sigs2)
0183 + x * ct * sin / sigs2
0184 + 2 * x * cos * cosks * sinkct * sinCR / k / sigmax2
0185
0186
0187 ) /
0188 (1.0 + (z * z) / (betx * betx)) / std::sqrt(1.0 + (z * z) / (bets * bets));
0189
0190 } else {
0191 const double norm = 2.0 / (NormFactorGauss4 * sigs2 * factorRMSgauss4);
0192 const double sigs4 = sigs2 * sigs2 * factorRMSgauss4 * factorRMSgauss4;
0193 const double cosks = std::cos(k * z);
0194 const double sinct = std::sin(k * ct);
0195 result =
0196 norm *
0197 std::exp(-ct * ct * ct * ct / sigs4 - z * z * z * z * cos2 * cos2 / sigs4 - 6 * ct * ct * z * z * cos2 / sigs4 -
0198 sin2 / (4 * k2 * sigmax2) *
0199 (2 + 4 * k2 * z * z - std::cos(2 * k * (z - ct)) - std::cos(2 * k * (z + ct)) -
0200 8 * k * s * std::cos(k * ct) * std::sin(k * z) - 4 * cosks * cosks * sinct * sinct)) /
0201 std::sqrt(1 + z * z / (betx * betx)) / std::sqrt(1 + z * z / (bets * bets));
0202 }
0203
0204 return result;
0205 }