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/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;            // eV
0021   constexpr double gamma34 = 1.22541670246517764513;  // Gamma(3/4)
0022   constexpr double gamma14 = 3.62560990822190831193;  // Gamma(1/4)
0023   constexpr double gamma54 = 0.90640247705547798267;  // Gamma(5/4)
0024   constexpr double sqrt2 = 1.41421356237;
0025   constexpr double sqrt2to5 = 5.65685424949;
0026   constexpr double two_pi = 2.0 * M_PI;
0027 }  // namespace
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   //---convert to mm
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   //---c in m/s --- remember t is already in meters
0128   constexpr double c = 2.99792458e+8;  // m/s
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;  // # Factor to take rms sigma as input of the supergaussian
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                               //-4*cosks*cosks * sinkct*sinkct * sinCR2 // comes from integral over x
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)               // contribution from x integrand
0183                       + x * ct * sin / sigs2                                // contribution from x integrand
0184                       + 2 * x * cos * cosks * sinkct * sinCR / k / sigmax2  // contribution from x integrand
0185                       //+(2*ct/k)*np.cos(k*s)*np.sin(k*ct) *(sin*sinCR)/(sigs2*cos)  # small term
0186                       //+ct**2*(sin2/sigs4)/(cos2/sigmax2)                            # small term
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 }