RivetInterface/plugins/MergedGenParticleProducer.cc

0001 #include <memory>
0002
0003 #include "MergedGenParticleProducer.hh"
0004
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008 #include "DataFormats/Common/interface/Handle.h"
0009 #include "DataFormats/Common/interface/View.h"
0010
0011 #include "HepPDT/ParticleID.hh"
0012
0013 MergedGenParticleProducer::MergedGenParticleProducer(const edm::ParameterSet& config) {
0014   input_pruned_ = consumes<edm::View<reco::GenParticle>>(config.getParameter<edm::InputTag>("inputPruned"));
0015   input_packed_ = consumes<edm::View<pat::PackedGenParticle>>(config.getParameter<edm::InputTag>("inputPacked"));
0016
0017   produces<reco::GenParticleCollection>();
0018 }
0019
0020 void MergedGenParticleProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
0021   // Need a ref to the product now for creating the mother/daughter refs
0022   auto ref = event.getRefBeforePut<reco::GenParticleCollection>();
0023
0024   // Get the input collections
0025   edm::Handle<edm::View<reco::GenParticle>> pruned_handle;
0026   event.getByToken(input_pruned_, pruned_handle);
0027
0028   edm::Handle<edm::View<pat::PackedGenParticle>> packed_handle;
0029   event.getByToken(input_packed_, packed_handle);
0030
0031   // First determine which packed particles are also still in the pruned collection
0032   // so that we can skip them later
0033   std::map<pat::PackedGenParticle const*, reco::GenParticle const*> st1_dup_map;
0034
0035   // Also map pointers in the original pruned collection to their index in the vector.
0036   // This index will be the same in the merged collection.
0037   std::map<reco::Candidate const*, std::size_t> pruned_idx_map;
0038
0039   unsigned int nLeptonsFromPrunedPhoton = 0;
0040
0041   for (unsigned int i = 0, idx = 0; i < pruned_handle->size(); ++i) {
0042     reco::GenParticle const& src = pruned_handle->at(i);
0043     pruned_idx_map[&src] = idx;
0044     ++idx;
0045
0046     // check for electrons+muons from pruned photons
0047     if (isLeptonFromPrunedPhoton(src)) {
0048       ++nLeptonsFromPrunedPhoton;
0049       ++idx;
0050     }
0051
0052     if (src.status() != 1)
0053       continue;
0054
0055     // Convert the pruned GenParticle into a PackedGenParticle then do an exact
0056     // floating point comparison of the pt, eta and phi
0057     // This of course relies on the PackedGenParticle constructor being identical
0058     // between the CMSSW version this sample was produced with and the one we're
0059     // analysing with
0060     pat::PackedGenParticle pks(src, reco::GenParticleRef());
0061     unsigned found_matches = 0;
0062     for (unsigned j = 0; j < packed_handle->size(); ++j) {
0063       pat::PackedGenParticle const& pk = packed_handle->at(j);
0064       if (pks.pdgId() != pk.pdgId() or pks.p4() != pk.p4())
0065         continue;
0066       ++found_matches;
0067       st1_dup_map[&pk] = &src;
0068     }
0069     if (found_matches > 1) {
0070       edm::LogWarning("MergedGenParticleProducer") << "Found multiple packed matches for: " << i << "\t" << src.pdgId()
0071                                                    << "\t" << src.pt() << "\t" << src.y() << "\n";
0072     } else if (found_matches == 0 && std::abs(src.y()) < 6.0) {
0073       edm::LogWarning("MergedGenParticleProducer")
0074           << "unmatched status 1: " << i << "\t" << src.pdgId() << "\t" << src.pt() << "\t" << src.y() << "\n";
0075     }
0076   }
0077
0078   // check for photons from pruned (light) hadrons
0079   unsigned int nPhotonsFromPrunedHadron = 0;
0080   for (unsigned int j = 0; j < packed_handle->size(); ++j) {
0081     pat::PackedGenParticle const& pk = packed_handle->at(j);
0082     if (isPhotonFromPrunedHadron(pk))
0083       ++nPhotonsFromPrunedHadron;
0084   }
0085
0086   // At this point we know what the size of the merged GenParticle will be so we can create it
0087   const unsigned int n = pruned_handle->size() + (packed_handle->size() - st1_dup_map.size()) +
0088                          nPhotonsFromPrunedHadron + nLeptonsFromPrunedPhoton;
0089   auto cands = std::make_unique<reco::GenParticleCollection>(n);
0090
0091   // First copy in all the pruned candidates
0092   unsigned idx = 0;
0093   for (unsigned i = 0; i < pruned_handle->size(); ++i) {
0094     reco::GenParticle const& old_cand = pruned_handle->at(i);
0095     reco::GenParticle& new_cand = cands->at(idx);
0096     new_cand = reco::GenParticle(pruned_handle->at(i));
0097     // Update the mother and daughter refs to this new merged collection
0098     new_cand.resetMothers(ref.id());
0099     new_cand.resetDaughters(ref.id());
0100     // Insert dummy photon mothers for orphaned electrons+muons
0101     if (isLeptonFromPrunedPhoton(old_cand)) {
0102       ++idx;
0103       reco::GenParticle& dummy_mother = cands->at(idx);
0104       dummy_mother = reco::GenParticle(0, old_cand.p4() * 1e-20, old_cand.vertex(), 22, 2, true);
0105       dummy_mother.addDaughter(reco::GenParticleRef(ref, idx - 1));
0106       new_cand.clearMothers();
0107       new_cand.addMother(reco::GenParticleRef(ref, idx));
0108       // now link the dummy to the packed candidate mothers
0109       for (unsigned m = 0; m < old_cand.numberOfMothers(); ++m) {
0110         dummy_mother.addMother(reco::GenParticleRef(ref, pruned_idx_map.at(old_cand.mother(m))));
0111         unsigned int midx = pruned_idx_map.at(old_cand.mother(m));
0112         if (midx < idx) {  // update existing mother to point to dummy
0113           reco::GenParticle& mother = cands->at(midx);
0114           mother.addDaughter(reco::GenParticleRef(ref, idx));
0115         } else {
0116           edm::LogWarning("MergedGenParticleProducer")
0117               << "Cannot assign to dummy photon with index " << idx
0118               << " as daughter to unprocessed particle with index " << idx << "\n";
0119         }
0120       }
0121     } else {
0122       for (unsigned m = 0; m < old_cand.numberOfMothers(); ++m) {
0123         new_cand.addMother(reco::GenParticleRef(ref, pruned_idx_map.at(old_cand.mother(m))));
0124       }
0125     }
0126     for (unsigned d = 0; d < old_cand.numberOfDaughters(); ++d) {
0127       new_cand.addDaughter(reco::GenParticleRef(ref, pruned_idx_map.at(old_cand.daughter(d))));
0128     }
0129     ++idx;
0130   }
0131
0132   // Now copy in the packed candidates that are not already in the pruned
0133   for (unsigned i = 0; i < packed_handle->size(); ++i) {
0134     pat::PackedGenParticle const& pk = packed_handle->at(i);
0135     if (st1_dup_map.count(&pk))
0136       continue;
0137     reco::GenParticle& new_cand = cands->at(idx);
0138     new_cand = reco::GenParticle(pk.charge(), pk.p4(), pk.vertex(), pk.pdgId(), 1, true);
0139
0140     // Insert dummy pi0 mothers for orphaned photons
0141     if (isPhotonFromPrunedHadron(pk)) {
0142       ++idx;
0143       reco::GenParticle& dummy_mother = cands->at(idx);
0144       dummy_mother = reco::GenParticle(0, pk.p4(), pk.vertex(), 111, 2, true);
0145       for (unsigned m = 0; m < pk.numberOfMothers(); ++m) {
0146         new_cand.addMother(reco::GenParticleRef(ref, idx));
0147         // Since the packed candidates drop the vertex position we'll take this from the mother
0148         if (m == 0) {
0149           dummy_mother.setP4(pk.mother(m)->p4());
0150           dummy_mother.setVertex(pk.mother(m)->vertex());
0151           new_cand.setVertex(pk.mother(m)->vertex());
0152         }
0153         // Should then add the GenParticle as a daughter of its dummy mother
0154         dummy_mother.addDaughter(reco::GenParticleRef(ref, idx - 1));
0155       }
0156     }
0157     // Connect to mother from pruned particles
0158     reco::GenParticle& daughter = cands->at(idx);
0159     for (unsigned m = 0; m < pk.numberOfMothers(); ++m) {
0160       daughter.addMother(reco::GenParticleRef(ref, pruned_idx_map.at(pk.mother(m))));
0161       // Since the packed candidates drop the vertex position we'll take this from the mother
0162       if (m == 0) {
0163         daughter.setVertex(pk.mother(m)->vertex());
0164       }
0165       // Should then add this GenParticle as a daughter of its mother
0166       cands->at(pruned_idx_map.at(pk.mother(m))).addDaughter(reco::GenParticleRef(ref, idx));
0167     }
0168     ++idx;
0169   }
0170
0171   event.put(std::move(cands));
0172 }
0173
0174 bool MergedGenParticleProducer::isPhotonFromPrunedHadron(const pat::PackedGenParticle& pk) const {
0175   if (pk.pdgId() == 22 and pk.statusFlags().isDirectHadronDecayProduct()) {
0176     // no mother
0177     if (pk.numberOfMothers() == 0)
0178       return true;
0179     // miniaod mother not compatible with the status flag
0180     HepPDT::ParticleID motherid(pk.mother(0)->pdgId());
0181     if (not(motherid.isHadron() and pk.mother(0)->status() == 2))
0182       return true;
0183   }
0184   return false;
0185 }
0186
0187 bool MergedGenParticleProducer::isLeptonFromPrunedPhoton(const reco::GenParticle& pk) const {
0188   if ((abs(pk.pdgId()) == 11 or abs(pk.pdgId()) == 13) and
0189       not(pk.statusFlags().fromHardProcess() or pk.statusFlags().isDirectTauDecayProduct())) {
0190     // this is probably not a prompt lepton but from pair production via a pruned photon
0191     if (pk.numberOfMothers() > 0 and pk.mother(0)->pdgId() != 22) {
0192       return true;
0193     }
0194   }
0195   return false;
0196 }
0197
0198 #include "FWCore/Framework/interface/MakerMacros.h"
0199 DEFINE_FWK_MODULE(MergedGenParticleProducer);