File indexing completed on 2025-01-28 23:24:34
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 #include "IOMC/EventVertexGenerators/interface/BetafuncEvtVtxGenerator.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023
0024 #include <CLHEP/Random/RandGaussQ.h>
0025 #include <CLHEP/Units/SystemOfUnits.h>
0026 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0027
0028 using CLHEP::cm;
0029 using CLHEP::ns;
0030 using CLHEP::radian;
0031
0032 BetafuncEvtVtxGenerator::BetafuncEvtVtxGenerator(const edm::ParameterSet& p) : BaseEvtVtxGenerator(p), boost_(4, 4) {
0033 readDB_ = p.getParameter<bool>("readDB");
0034 if (!readDB_) {
0035 fX0 = p.getParameter<double>("X0") * cm;
0036 fY0 = p.getParameter<double>("Y0") * cm;
0037 fZ0 = p.getParameter<double>("Z0") * cm;
0038 fSigmaZ = p.getParameter<double>("SigmaZ") * cm;
0039 fbetastar = p.getParameter<double>("BetaStar") * cm;
0040 femittance = p.getParameter<double>("Emittance") * cm;
0041 fTimeOffset = p.getParameter<double>("TimeOffset") * ns * c_light;
0042
0043 setBoost(p.getParameter<double>("Alpha") * radian, p.getParameter<double>("Phi") * radian);
0044 if (fSigmaZ <= 0) {
0045 throw cms::Exception("Configuration") << "Error in BetafuncEvtVtxGenerator: "
0046 << "Illegal resolution in Z (SigmaZ is negative)";
0047 }
0048 }
0049 if (readDB_) {
0050
0051
0052 beamToken_ = esConsumes<SimBeamSpotObjects, SimBeamSpotObjectsRcd, edm::Transition::BeginLuminosityBlock>();
0053 }
0054 }
0055
0056 void BetafuncEvtVtxGenerator::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const& iEventSetup) {
0057 update(iEventSetup);
0058 }
0059
0060 void BetafuncEvtVtxGenerator::update(const edm::EventSetup& iEventSetup) {
0061 if (readDB_ && parameterWatcher_.check(iEventSetup)) {
0062 edm::ESHandle<SimBeamSpotObjects> beamhandle = iEventSetup.getHandle(beamToken_);
0063 if (!beamhandle->isGaussian()) {
0064 fX0 = beamhandle->x() * cm;
0065 fY0 = beamhandle->y() * cm;
0066 fZ0 = beamhandle->z() * cm;
0067 fSigmaZ = beamhandle->sigmaZ() * cm;
0068 fTimeOffset = beamhandle->timeOffset() * ns * c_light;
0069 fbetastar = beamhandle->betaStar() * cm;
0070 femittance = beamhandle->emittance() * cm;
0071 setBoost(beamhandle->alpha() * radian, beamhandle->phi() * radian);
0072 } else {
0073 throw cms::Exception("Configuration")
0074 << "Error in BetafuncEvtVtxGenerator::update: The provided SimBeamSpotObjects is Gaussian.\n"
0075 << "Please check the configuration and ensure that the beam spot parameters are appropriate for a Betafunc "
0076 "distribution.";
0077 }
0078 }
0079 }
0080
0081 ROOT::Math::XYZTVector BetafuncEvtVtxGenerator::vertexShift(CLHEP::HepRandomEngine* engine) const {
0082 double X, Y, Z;
0083
0084 double tmp_sigz = CLHEP::RandGaussQ::shoot(engine, 0., fSigmaZ);
0085 Z = tmp_sigz + fZ0;
0086
0087 double tmp_sigx = BetaFunction(Z, fZ0);
0088
0089 tmp_sigx /= sqrt(2.0);
0090 X = CLHEP::RandGaussQ::shoot(engine, 0., tmp_sigx) + fX0;
0091
0092 double tmp_sigy = BetaFunction(Z, fZ0);
0093
0094 tmp_sigy /= sqrt(2.0);
0095 Y = CLHEP::RandGaussQ::shoot(engine, 0., tmp_sigy) + fY0;
0096
0097 double tmp_sigt = CLHEP::RandGaussQ::shoot(engine, 0., fSigmaZ);
0098 double T = tmp_sigt + fTimeOffset;
0099
0100 return ROOT::Math::XYZTVector(X, Y, Z, T);
0101 }
0102
0103 double BetafuncEvtVtxGenerator::BetaFunction(double z, double z0) const {
0104 return sqrt(femittance * (fbetastar + (((z - z0) * (z - z0)) / fbetastar)));
0105 }
0106
0107 void BetafuncEvtVtxGenerator::setBoost(double alpha, double phi) {
0108
0109
0110 TMatrixD tmpboost(4, 4);
0111
0112
0113
0114
0115
0116
0117
0118 tmpboost(0, 0) = 1. / cos(phi);
0119 tmpboost(0, 1) = -cos(alpha) * sin(phi);
0120 tmpboost(0, 2) = -tan(phi) * sin(phi);
0121 tmpboost(0, 3) = -sin(alpha) * sin(phi);
0122 tmpboost(1, 0) = -cos(alpha) * tan(phi);
0123 tmpboost(1, 1) = 1.;
0124 tmpboost(1, 2) = cos(alpha) * tan(phi);
0125 tmpboost(1, 3) = 0.;
0126 tmpboost(2, 0) = 0.;
0127 tmpboost(2, 1) = -cos(alpha) * sin(phi);
0128 tmpboost(2, 2) = cos(phi);
0129 tmpboost(2, 3) = -sin(alpha) * sin(phi);
0130 tmpboost(3, 0) = -sin(alpha) * tan(phi);
0131 tmpboost(3, 1) = 0.;
0132 tmpboost(3, 2) = sin(alpha) * tan(phi);
0133 tmpboost(3, 3) = 1.;
0134
0135 tmpboost.Invert();
0136 boost_ = tmpboost;
0137
0138 }
0139
0140 void BetafuncEvtVtxGenerator::sigmaZ(double s) {
0141 if (s >= 0) {
0142 fSigmaZ = s;
0143 } else {
0144 throw cms::Exception("LogicError") << "Error in BetafuncEvtVtxGenerator::sigmaZ: "
0145 << "Illegal resolution in Z (negative)";
0146 }
0147 }
0148
0149 TMatrixD const* BetafuncEvtVtxGenerator::GetInvLorentzBoost() const { return &boost_; }
0150
0151 void BetafuncEvtVtxGenerator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0152 edm::ParameterSetDescription desc;
0153 desc.add<double>("X0", 0.0)->setComment("in cm");
0154 desc.add<double>("Y0", 0.0)->setComment("in cm");
0155 desc.add<double>("Z0", 0.0)->setComment("in cm");
0156 desc.add<double>("SigmaZ", 0.0)->setComment("in cm");
0157 desc.add<double>("BetaStar", 0.0)->setComment("in cm");
0158 desc.add<double>("Emittance", 0.0)->setComment("in cm");
0159 desc.add<double>("Alpha", 0.0)->setComment("in radians");
0160 desc.add<double>("Phi", 0.0)->setComment("in radians");
0161 desc.add<double>("TimeOffset", 0.0)->setComment("in ns");
0162 desc.add<edm::InputTag>("src");
0163 desc.add<bool>("readDB");
0164 descriptions.add("BetafuncEvtVtxGenerator", desc);
0165 }