Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:21:07

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);