GenParticleDecaySelector.cc

Go to the documentation of this file.

/* \class GenParticleDecaySelector
 *
 * \author Luca Lista, INFN
 *
 *
 */

#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "SimGeneral/HepPDTRecord/interface/PdtEntry.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"

class GenParticleDecaySelector : public edm::EDProducer {
public:
  GenParticleDecaySelector(const edm::ParameterSet&);

private:
  void produce(edm::Event& e, const edm::EventSetup&) override;
  bool firstEvent_;
  edm::EDGetTokenT<reco::GenParticleCollection> srcToken_;
  PdtEntry particle_;
  int status_;
  std::pair<reco::GenParticleRef, reco::GenParticle*> add(reco::GenParticleCollection&,
                                                          const reco::GenParticle&,
                                                          reco::GenParticleRefProd);
};

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/Event.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"

using namespace edm;
using namespace reco;
using namespace std;

GenParticleDecaySelector::GenParticleDecaySelector(const edm::ParameterSet& cfg)
    : firstEvent_(true),
      srcToken_(consumes<GenParticleCollection>(cfg.getParameter<InputTag>("src"))),
      particle_(cfg.getParameter<PdtEntry>("particle")),
      status_(cfg.getParameter<int>("status")) {
  produces<GenParticleCollection>();
}

void GenParticleDecaySelector::produce(edm::Event& evt, const edm::EventSetup& es) {
  if (firstEvent_) {
    particle_.setup(es);
    firstEvent_ = false;
  }

  Handle<GenParticleCollection> genParticles;
  evt.getByToken(srcToken_, genParticles);
  auto decay = std::make_unique<GenParticleCollection>();
  const GenParticleRefProd ref = evt.getRefBeforePut<GenParticleCollection>();
  for (GenParticleCollection::const_iterator g = genParticles->begin(); g != genParticles->end(); ++g)
    if (g->pdgId() == particle_.pdgId() && g->status() == status_)
      add(*decay, *g, ref);
  evt.put(std::move(decay));
}

pair<GenParticleRef, GenParticle*> GenParticleDecaySelector::add(GenParticleCollection& decay,
                                                                 const GenParticle& p,
                                                                 GenParticleRefProd ref) {
  size_t idx = decay.size();
  GenParticleRef r(ref, idx);
  decay.resize(idx + 1);
  const LeafCandidate& part = p;
  GenParticle g(part);
  size_t n = p.numberOfDaughters();
  for (size_t i = 0; i < n; ++i) {
    pair<GenParticleRef, GenParticle*> d = add(decay, *p.daughterRef(i), ref);
    d.second->addMother(r);
    g.addDaughter(d.first);
  }
  GenParticle& gp = decay[idx] = g;
  return make_pair(r, &gp);
}

#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(GenParticleDecaySelector);