Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-31 22:54:36

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   // Fix by Markus
0079   // check for photons from pruned (light) hadrons
0080   unsigned int nPhotonsFromPrunedHadron = 0;
0081   for (unsigned int j = 0; j < packed_handle->size(); ++j) {
0082     pat::PackedGenParticle const& pk = packed_handle->at(j);
0083     if (isPhotonFromPrunedHadron(pk))
0084       ++nPhotonsFromPrunedHadron;
0085   }
0086 
0087   // At this point we know what the size of the merged GenParticle will be so we can create it
0088   const unsigned int n = pruned_handle->size() + (packed_handle->size() - st1_dup_map.size()) +
0089                          nPhotonsFromPrunedHadron + nLeptonsFromPrunedPhoton;
0090   auto cands = std::make_unique<reco::GenParticleCollection>(n);
0091 
0092   // First copy in all the pruned candidates
0093   unsigned idx = 0;
0094   for (unsigned i = 0; i < pruned_handle->size(); ++i) {
0095     reco::GenParticle const& old_cand = pruned_handle->at(i);
0096     reco::GenParticle& new_cand = cands->at(idx);
0097     new_cand = reco::GenParticle(pruned_handle->at(i));
0098     // Update the mother and daughter refs to this new merged collection
0099     new_cand.resetMothers(ref.id());
0100     new_cand.resetDaughters(ref.id());
0101     // Insert dummy photon mothers for orphaned electrons+muons
0102     if (isLeptonFromPrunedPhoton(old_cand)) {
0103       ++idx;
0104       reco::GenParticle& dummy_mother = cands->at(idx);
0105       dummy_mother = reco::GenParticle(0, old_cand.p4(), old_cand.vertex(), 22, 2, true);
0106       for (unsigned m = 0; m < old_cand.numberOfMothers(); ++m) {
0107         new_cand.addMother(reco::GenParticleRef(ref, idx));
0108         // Since the packed candidates drop the vertex position we'll take this from the mother
0109         if (m == 0) {
0110           dummy_mother.setP4(old_cand.mother(0)->p4());
0111           dummy_mother.setVertex(old_cand.mother(0)->vertex());
0112           new_cand.setVertex(old_cand.mother(0)->vertex());
0113         }
0114         // Should then add the GenParticle as a daughter of its dummy mother
0115         dummy_mother.addDaughter(reco::GenParticleRef(ref, idx - 1));
0116         for (unsigned m = 0; m < old_cand.numberOfMothers(); ++m) {
0117           dummy_mother.addMother(reco::GenParticleRef(ref, pruned_idx_map.at(old_cand.mother(m))));
0118         }
0119       }
0120     } else {
0121       for (unsigned m = 0; m < old_cand.numberOfMothers(); ++m) {
0122         new_cand.addMother(reco::GenParticleRef(ref, pruned_idx_map.at(old_cand.mother(m))));
0123       }
0124     }
0125     for (unsigned d = 0; d < old_cand.numberOfDaughters(); ++d) {
0126       new_cand.addDaughter(reco::GenParticleRef(ref, pruned_idx_map.at(old_cand.daughter(d))));
0127     }
0128     ++idx;
0129   }
0130 
0131   // Now copy in the packed candidates that are not already in the pruned
0132   for (unsigned i = 0; i < packed_handle->size(); ++i) {
0133     pat::PackedGenParticle const& pk = packed_handle->at(i);
0134     if (st1_dup_map.count(&pk))
0135       continue;
0136     reco::GenParticle& new_cand = cands->at(idx);
0137     new_cand = reco::GenParticle(pk.charge(), pk.p4(), pk.vertex(), pk.pdgId(), 1, true);
0138 
0139     // Insert dummy pi0 mothers for orphaned photons
0140     if (isPhotonFromPrunedHadron(pk)) {
0141       ++idx;
0142       reco::GenParticle& dummy_mother = cands->at(idx);
0143       dummy_mother = reco::GenParticle(0, pk.p4(), pk.vertex(), 111, 2, true);
0144       for (unsigned m = 0; m < pk.numberOfMothers(); ++m) {
0145         new_cand.addMother(reco::GenParticleRef(ref, idx));
0146         // Since the packed candidates drop the vertex position we'll take this from the mother
0147         if (m == 0) {
0148           dummy_mother.setP4(pk.mother(m)->p4());
0149           dummy_mother.setVertex(pk.mother(m)->vertex());
0150           new_cand.setVertex(pk.mother(m)->vertex());
0151         }
0152         // Should then add the GenParticle as a daughter of its dummy mother
0153         dummy_mother.addDaughter(reco::GenParticleRef(ref, idx - 1));
0154       }
0155     }
0156     // Connect to mother from pruned particles
0157     reco::GenParticle& daughter = cands->at(idx);
0158     for (unsigned m = 0; m < pk.numberOfMothers(); ++m) {
0159       daughter.addMother(reco::GenParticleRef(ref, pruned_idx_map.at(pk.mother(m))));
0160       // Since the packed candidates drop the vertex position we'll take this from the mother
0161       if (m == 0) {
0162         daughter.setVertex(pk.mother(m)->vertex());
0163       }
0164       // Should then add this GenParticle as a daughter of its mother
0165       cands->at(pruned_idx_map.at(pk.mother(m))).addDaughter(reco::GenParticleRef(ref, idx));
0166     }
0167     ++idx;
0168   }
0169 
0170   event.put(std::move(cands));
0171 }
0172 
0173 bool MergedGenParticleProducer::isPhotonFromPrunedHadron(const pat::PackedGenParticle& pk) const {
0174   if (pk.pdgId() == 22 and pk.statusFlags().isDirectHadronDecayProduct()) {
0175     // no mother
0176     if (pk.numberOfMothers() == 0)
0177       return true;
0178     // miniaod mother not compatible with the status flag
0179     HepPDT::ParticleID motherid(pk.mother(0)->pdgId());
0180     if (not(motherid.isHadron() and pk.mother(0)->status() == 2))
0181       return true;
0182   }
0183   return false;
0184 }
0185 
0186 bool MergedGenParticleProducer::isLeptonFromPrunedPhoton(const reco::GenParticle& pk) const {
0187   if ((abs(pk.pdgId()) == 11 or abs(pk.pdgId()) == 13) and
0188       not(pk.statusFlags().fromHardProcess() or pk.statusFlags().isDirectTauDecayProduct())) {
0189     // this is probably not a prompt lepton but from pair production via a pruned photon
0190     if (pk.numberOfMothers() > 0 and pk.mother(0)->pdgId() != 22) {
0191       return true;
0192     }
0193   }
0194   return false;
0195 }
0196 
0197 #include "FWCore/Framework/interface/MakerMacros.h"
0198 DEFINE_FWK_MODULE(MergedGenParticleProducer);