00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "FWCore/Framework/interface/EDProducer.h"
00010 #include "FWCore/ParameterSet/interface/InputTag.h"
00011 #include "SimGeneral/HepPDTRecord/interface/PdtEntry.h"
00012 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
00013
00014 class GenParticleDecaySelector : public edm::EDProducer {
00015 public:
00017 GenParticleDecaySelector(const edm::ParameterSet&);
00018 private:
00020 void produce(edm::Event& e, const edm::EventSetup&);
00022 void beginJob(const edm::EventSetup&);
00024 edm::InputTag src_;
00026 PdtEntry particle_;
00028 int status_;
00030 std::pair<reco::GenParticleRef, reco::GenParticle*>
00031 add(reco::GenParticleCollection&, const reco::GenParticle &,
00032 reco::GenParticleRefProd);
00033 };
00034
00035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00036 #include "DataFormats/Common/interface/Handle.h"
00037 #include "FWCore/Framework/interface/Event.h"
00038 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00039
00040 using namespace edm;
00041 using namespace reco;
00042 using namespace std;
00043
00044 GenParticleDecaySelector::GenParticleDecaySelector(const edm::ParameterSet& cfg) :
00045 src_(cfg.getParameter<InputTag>("src")),
00046 particle_(cfg.getParameter<PdtEntry>("particle")),
00047 status_(cfg.getParameter<int>("status")) {
00048 produces<GenParticleCollection>();
00049 }
00050
00051 void GenParticleDecaySelector::beginJob(const edm::EventSetup &es) {
00052 particle_.setup(es);
00053 }
00054
00055 void GenParticleDecaySelector::produce(edm::Event& evt, const edm::EventSetup&) {
00056 Handle<GenParticleCollection> genParticles;
00057 evt.getByLabel(src_, genParticles);
00058 auto_ptr<GenParticleCollection> decay(new GenParticleCollection);
00059 const GenParticleRefProd ref = evt.getRefBeforePut<GenParticleCollection>();
00060 for(GenParticleCollection::const_iterator g = genParticles->begin();
00061 g != genParticles->end(); ++g)
00062 if(g->pdgId() == particle_.pdgId() && g->status() == status_)
00063 add(*decay, *g, ref);
00064 evt.put(decay);
00065 }
00066
00067 pair<GenParticleRef, GenParticle*> GenParticleDecaySelector::add(GenParticleCollection & decay, const GenParticle & p,
00068 GenParticleRefProd ref) {
00069 size_t idx = decay.size();
00070 GenParticleRef r(ref, idx);
00071 decay.resize(idx+1);
00072 const Particle & part = p;
00073 GenParticle g(part);
00074 size_t n = p.numberOfDaughters();
00075 for(size_t i = 0; i < n; ++i) {
00076 pair<GenParticleRef, GenParticle*> d = add(decay, *p.daughterRef(i), ref);
00077 d.second->addMother(r);
00078 g.addDaughter(d.first);
00079 }
00080 GenParticle & gp = decay[idx] = g;
00081 return make_pair(r, &gp);
00082 }
00083
00084
00085 #include "FWCore/Framework/interface/MakerMacros.h"
00086
00087 DEFINE_FWK_MODULE(GenParticleDecaySelector);