EventVertexGenerators/src/HLLHCEvtVtxGenerator.cc

0001 // $Id: HLLHCEvtVtxGenerator_Fix.cc, v 1.0 2015/03/15 10:32:11 Exp $
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;            // eV
0019   constexpr double gamma34 = 1.22541670246517764513;  // Gamma(3/4)
0020   constexpr double gamma14 = 3.62560990822190831193;  // Gamma(1/4)
0021   constexpr double gamma54 = 0.90640247705547798267;  // Gamma(5/4)
0022   constexpr double sqrt2 = 1.41421356237;
0023   constexpr double sqrt2to5 = 5.65685424949;
0024   constexpr double two_pi = 2.0 * M_PI;
0025 }  // namespace
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     // Read configurable parameters
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     // Set parameters inferred from configurables
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     // NOTE: this is currently watching LS transitions, while it should watch Run transitions,
0080     // even though in reality there is no Run Dependent MC (yet) in CMS
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     // Read configurable parameters
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     // Set parameters inferred from configurables
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   //---convert to mm
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   //---c in m/s --- remember t is already in meters
0167   constexpr double c = 2.99792458e+8;  // m/s
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;  // # Factor to take rms sigma as input of the supergaussian
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                               //-4*cosks*cosks * sinkct*sinkct * sinCR2 // comes from integral over x
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)               // contribution from x integrand
0222                       + x * ct * sin / sigs2                                // contribution from x integrand
0223                       + 2 * x * cos * cosks * sinkct * sinCR / k / sigmax2  // contribution from x integrand
0224                       //+(2*ct/k)*np.cos(k*s)*np.sin(k*ct) *(sin*sinCR)/(sigs2*cos)  # small term
0225                       //+ct**2*(sin2/sigs4)/(cos2/sigmax2)                           # small term
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 }