EventVertexGenerators/src/BeamDivergenceVtxGenerator.cc

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  ****************************************************************************/
0005
0006 #include "FWCore/Framework/interface/stream/EDProducer.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0016 #include "FWCore/Utilities/interface/ESGetToken.h"
0017
0018 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0019
0020 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0021
0022 #include "CondFormats/PPSObjects/interface/CTPPSBeamParameters.h"
0023 #include "CondFormats/DataRecord/interface/CTPPSBeamParametersRcd.h"
0024
0025 #include <CLHEP/Random/RandGauss.h>
0026
0027 //----------------------------------------------------------------------------------------------------
0028
0029 class BeamDivergenceVtxGenerator : public edm::stream::EDProducer<> {
0030 public:
0031   explicit BeamDivergenceVtxGenerator(const edm::ParameterSet &);
0032   ~BeamDivergenceVtxGenerator() override = default;
0033
0034   static void fillDescriptions(edm::ConfigurationDescriptions &);
0035
0036   void produce(edm::Event &, const edm::EventSetup &) override;
0037
0038 private:
0039   edm::EDGetTokenT<edm::HepMCProduct> sourceToken_;
0040   edm::ESGetToken<CTPPSBeamParameters, CTPPSBeamParametersRcd> beamParametersToken_;
0041   std::vector<edm::EDGetTokenT<reco::GenParticleCollection>> genParticleTokens_;
0042
0043   bool simulateVertex_;
0044   bool simulateBeamDivergence_;
0045
0046   struct SmearingParameters {
0047     double vtx_x, vtx_y, vtx_z, vtx_t;          // cm
0048     double bd_x_45, bd_y_45, bd_x_56, bd_y_56;  // rad
0049   };
0050
0051   template <typename T>
0052   static HepMC::FourVector smearMomentum(const T &mom, const SmearingParameters &sp);
0053
0054   void applySmearingHepMC(const SmearingParameters &sp, HepMC::GenEvent *genEvt);
0055
0056   void addSmearedGenParticle(const reco::GenParticle &gp, const SmearingParameters &sp, HepMC::GenEvent *genEvt);
0057 };
0058
0059 //----------------------------------------------------------------------------------------------------
0060
0061 BeamDivergenceVtxGenerator::BeamDivergenceVtxGenerator(const edm::ParameterSet &iConfig)
0062     : beamParametersToken_(esConsumes<CTPPSBeamParameters, CTPPSBeamParametersRcd>()),
0063       simulateVertex_(iConfig.getParameter<bool>("simulateVertex")),
0064       simulateBeamDivergence_(iConfig.getParameter<bool>("simulateBeamDivergence")) {
0065   const edm::InputTag tagSrcHepMC = iConfig.getParameter<edm::InputTag>("src");
0066   if (!tagSrcHepMC.label().empty())
0067     sourceToken_ = consumes<edm::HepMCProduct>(tagSrcHepMC);
0068
0069   for (const auto &it : iConfig.getParameter<std::vector<edm::InputTag>>("srcGenParticle"))
0070     genParticleTokens_.push_back(consumes<reco::GenParticleCollection>(it));
0071
0072   edm::Service<edm::RandomNumberGenerator> rng;
0073   if (!rng.isAvailable())
0074     throw cms::Exception("Configuration")
0075         << "The BeamDivergenceVtxGenerator requires the RandomNumberGeneratorService\n"
0076            "which is not present in the configuration file. \n"
0077            "You must add the service\n"
0078            "in the configuration file or remove the modules that require it.";
0079
0080   produces<edm::HepMCProduct>();
0081 }
0082
0083 //----------------------------------------------------------------------------------------------------
0084
0085 void BeamDivergenceVtxGenerator::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0086   edm::ParameterSetDescription desc;
0087
0088   desc.add<edm::InputTag>("src", edm::InputTag("generator", "unsmeared"))
0089       ->setComment("input collection in HepMC format");
0090
0091   desc.add<std::vector<edm::InputTag>>("srcGenParticle", std::vector<edm::InputTag>())
0092       ->setComment("input collections in GenParticle format");
0093
0094   desc.add<bool>("simulateBeamDivergence", true)->setComment("account for the beam angular divergence?");
0095   desc.add<bool>("simulateVertex", true)->setComment("account for the vertex transverse smearing?");
0096
0097   descriptions.add("beamDivergenceVtxGenerator", desc);
0098 }
0099
0100 //----------------------------------------------------------------------------------------------------
0101
0102 void BeamDivergenceVtxGenerator::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0103   // get random engine
0104   edm::Service<edm::RandomNumberGenerator> rng;
0105   CLHEP::HepRandomEngine *rnd = &(rng->getEngine(iEvent.streamID()));
0106
0107   // get conditions
0108   edm::ESHandle<CTPPSBeamParameters> hBeamParameters = iSetup.getHandle(beamParametersToken_);
0109
0110   // get HepMC input (if given)
0111   HepMC::GenEvent *genEvt;
0112   if (sourceToken_.isUninitialized()) {
0113     genEvt = new HepMC::GenEvent();
0114   } else {
0115     edm::Handle<edm::HepMCProduct> hepUnsmearedMCEvt;
0116     iEvent.getByToken(sourceToken_, hepUnsmearedMCEvt);
0117
0118     genEvt = new HepMC::GenEvent(*hepUnsmearedMCEvt->GetEvent());
0119   }
0120
0121   // prepare output
0122   std::unique_ptr<edm::HepMCProduct> output(new edm::HepMCProduct(genEvt));
0123
0124   // generate smearing parameters
0125   SmearingParameters sp;
0126
0127   if (simulateVertex_) {
0128     // NB: the separtion between effective offsets in LHC sectors 45 and 56 cannot be applied, thus the values for 45 are used
0129     sp.vtx_x = hBeamParameters->getVtxOffsetX45() + CLHEP::RandGauss::shoot(rnd) * hBeamParameters->getVtxStddevX();
0130     sp.vtx_y = hBeamParameters->getVtxOffsetY45() + CLHEP::RandGauss::shoot(rnd) * hBeamParameters->getVtxStddevY();
0131     sp.vtx_z = hBeamParameters->getVtxOffsetZ45() + CLHEP::RandGauss::shoot(rnd) * hBeamParameters->getVtxStddevZ();
0132     sp.vtx_t = hBeamParameters->getVtxOffsetT45() + CLHEP::RandGauss::shoot(rnd) * hBeamParameters->getVtxStddevT();
0133   }
0134
0135   if (simulateBeamDivergence_) {
0136     sp.bd_x_45 = CLHEP::RandGauss::shoot(rnd) * hBeamParameters->getBeamDivergenceX45();
0137     sp.bd_x_56 = CLHEP::RandGauss::shoot(rnd) * hBeamParameters->getBeamDivergenceX56();
0138     sp.bd_y_45 = CLHEP::RandGauss::shoot(rnd) * hBeamParameters->getBeamDivergenceY45();
0139     sp.bd_y_56 = CLHEP::RandGauss::shoot(rnd) * hBeamParameters->getBeamDivergenceY56();
0140   }
0141
0142   // apply smearing
0143   applySmearingHepMC(sp, genEvt);
0144
0145   for (const auto &token : genParticleTokens_) {
0146     edm::Handle<reco::GenParticleCollection> hGPCollection;
0147     iEvent.getByToken(token, hGPCollection);
0148
0149     for (const auto &gp : *hGPCollection)
0150       addSmearedGenParticle(gp, sp, genEvt);
0151   }
0152
0153   // save output
0154   iEvent.put(std::move(output));
0155 }
0156
0157 //----------------------------------------------------------------------------------------------------
0158
0159 template <typename T>
0160 HepMC::FourVector BeamDivergenceVtxGenerator::smearMomentum(const T &mom, const SmearingParameters &sp) {
0161   // TODO: this is an oversimplified implemetation
0162   // the TOTEM smearing module should be taken as reference
0163
0164   double th_x = mom.x() / mom.z();
0165   double th_y = mom.y() / mom.z();
0166
0167   if (mom.z() > 0.) {
0168     th_x += sp.bd_x_45;
0169     th_y += sp.bd_y_45;
0170   } else {
0171     th_x += sp.bd_x_56;
0172     th_y += sp.bd_y_56;
0173   }
0174
0175   // calculate consistent p_z component
0176   const double sign = (mom.z() > 0.) ? 1. : -1.;
0177   const double p = sqrt(mom.x() * mom.x() + mom.y() * mom.y() + mom.z() * mom.z());
0178   const double p_z = sign * p / sqrt(1. + th_x * th_x + th_y * th_y);
0179
0180   return HepMC::FourVector(p_z * th_x, p_z * th_y, p_z, mom.e());
0181 }
0182
0183 //----------------------------------------------------------------------------------------------------
0184
0185 void BeamDivergenceVtxGenerator::applySmearingHepMC(const SmearingParameters &sp, HepMC::GenEvent *genEvt) {
0186   if (simulateVertex_) {
0187     for (HepMC::GenEvent::vertex_iterator vit = genEvt->vertices_begin(); vit != genEvt->vertices_end(); ++vit) {
0188       const auto &pos = (*vit)->position();
0189       (*vit)->set_position(HepMC::FourVector(pos.x() + sp.vtx_x * 1E1,  // conversion: cm to mm
0190                                              pos.y() + sp.vtx_y * 1E1,
0191                                              pos.z() + sp.vtx_z * 1E1,
0192                                              pos.t() + sp.vtx_t * 1E1));
0193     }
0194   }
0195
0196   if (simulateBeamDivergence_) {
0197     for (HepMC::GenEvent::particle_iterator part = genEvt->particles_begin(); part != genEvt->particles_end(); ++part)
0198       (*part)->set_momentum(smearMomentum((*part)->momentum(), sp));
0199   }
0200 }
0201
0202 //----------------------------------------------------------------------------------------------------
0203
0204 void BeamDivergenceVtxGenerator::addSmearedGenParticle(const reco::GenParticle &gp,
0205                                                        const SmearingParameters &sp,
0206                                                        HepMC::GenEvent *genEvt) {
0207   // add vertex of the particle
0208   HepMC::GenVertex *vtx = new HepMC::GenVertex(HepMC::FourVector(
0209       (gp.vx() + sp.vtx_x) * 1E1,  // conversion: cm to mm
0210       (gp.vy() + sp.vtx_y) * 1E1,
0211       (gp.vz() + sp.vtx_z) * 1E1,
0212       (/*gp.vt()*/ +sp.vtx_t) * 1E1));  // TODO: GenParticle doesn't seem to have time component of the vertex
0213   genEvt->add_vertex(vtx);
0214
0215   // add the particle itself
0216   HepMC::GenParticle *particle = new HepMC::GenParticle(smearMomentum(gp.p4(), sp), gp.pdgId(), gp.status());
0217   vtx->add_particle_out(particle);
0218 }
0219
0220 //----------------------------------------------------------------------------------------------------
0221
0222 DEFINE_FWK_MODULE(BeamDivergenceVtxGenerator);