-
Notifications
You must be signed in to change notification settings - Fork 4.2k
/
MergedGenParticleProducer.cc
146 lines (126 loc) · 6.48 KB
/
MergedGenParticleProducer.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
#include "MergedGenParticleProducer.hh"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Common/interface/View.h"
#include "HepPDT/ParticleID.hh"
MergedGenParticleProducer::MergedGenParticleProducer(const edm::ParameterSet& config)
{
input_pruned_ = consumes<edm::View<reco::GenParticle>>(config.getParameter<edm::InputTag>("inputPruned"));
input_packed_ = consumes<edm::View<pat::PackedGenParticle>>(config.getParameter<edm::InputTag>("inputPacked"));
produces<reco::GenParticleCollection>();
}
void MergedGenParticleProducer::produce(edm::Event& event, const edm::EventSetup& setup)
{
// Need a ref to the product now for creating the mother/daughter refs
auto ref = event.getRefBeforePut<reco::GenParticleCollection>();
// Get the input collections
edm::Handle<edm::View<reco::GenParticle> > pruned_handle;
event.getByToken(input_pruned_, pruned_handle);
edm::Handle<edm::View<pat::PackedGenParticle> > packed_handle;
event.getByToken(input_packed_, packed_handle);
// First determine which packed particles are also still in the pruned collection
// so that we can skip them later
std::map<pat::PackedGenParticle const*, reco::GenParticle const*> st1_dup_map;
// Also map pointers in the original pruned collection to their index in the vector.
// This index will be the same in the merged collection.
std::map<reco::Candidate const*, std::size_t> pruned_idx_map;
for (unsigned int i = 0; i < pruned_handle->size(); ++i) {
reco::GenParticle const& src = pruned_handle->at(i);
pruned_idx_map[&src] = i;
if (src.status() != 1) continue;
// Convert the pruned GenParticle into a PackedGenParticle then do an exact
// floating point comparison of the pt, eta and phi
// This of course relies on the PackedGenParticle constructor being identical
// between the CMSSW version this sample was produced with and the one we're
// analysing with
pat::PackedGenParticle pks(src, reco::GenParticleRef());
unsigned found_matches = 0;
for (unsigned j = 0; j < packed_handle->size(); ++j) {
pat::PackedGenParticle const& pk = packed_handle->at(j);
if ( pks.pdgId() != pk.pdgId() or pks.p4() != pk.p4() ) continue;
++found_matches;
st1_dup_map[&pk] = &src;
}
if (found_matches > 1) {
edm::LogWarning("MergedGenParticleProducer") << "Found multiple packed matches for: " << i << "\t" << src.pdgId() << "\t" << src.pt() << "\t" << src.y() << "\n";
}
else if (found_matches == 0 && std::abs(src.y()) < 6.0) {
edm::LogWarning("MergedGenParticleProducer") << "unmatched status 1: " << i << "\t" << src.pdgId() << "\t" << src.pt() << "\t" << src.y() << "\n";
}
}
// Fix by Markus
// check for photons from pruned (light) hadrons
unsigned int nPhotonsFromPrunedHadron = 0;
for (unsigned int j = 0; j < packed_handle->size(); ++j) {
pat::PackedGenParticle const& pk = packed_handle->at(j);
if (isPhotonFromPrunedHadron(pk)) ++nPhotonsFromPrunedHadron;
}
// At this point we know what the size of the merged GenParticle will be so we can create it
const unsigned int n = pruned_handle->size() + (packed_handle->size() - st1_dup_map.size()) + nPhotonsFromPrunedHadron;
auto cands = std::unique_ptr<reco::GenParticleCollection>(new reco::GenParticleCollection(n));
// First copy in all the pruned candidates
for (unsigned i = 0; i < pruned_handle->size(); ++i) {
reco::GenParticle const& old_cand = pruned_handle->at(i);
reco::GenParticle & new_cand = cands->at(i);
new_cand = reco::GenParticle(pruned_handle->at(i));
// Update the mother and daughter refs to this new merged collection
new_cand.resetMothers(ref.id());
new_cand.resetDaughters(ref.id());
for (unsigned m = 0; m < old_cand.numberOfMothers(); ++m) {
new_cand.addMother(reco::GenParticleRef(ref, pruned_idx_map.at(old_cand.mother(m))));
}
for (unsigned d = 0; d < old_cand.numberOfDaughters(); ++d) {
new_cand.addDaughter(reco::GenParticleRef(ref, pruned_idx_map.at(old_cand.daughter(d))));
}
}
// Now copy in the packed candidates that are not already in the pruned
for (unsigned i = 0, idx = pruned_handle->size(); i < packed_handle->size(); ++i) {
pat::PackedGenParticle const& pk = packed_handle->at(i);
if (st1_dup_map.count(&pk)) continue;
reco::GenParticle & new_cand = cands->at(idx);
new_cand = reco::GenParticle(pk.charge(), pk.p4(), pk.vertex(), pk.pdgId(), 1, true);
// Insert dummy pi0 mothers for orphaned photons
if (isPhotonFromPrunedHadron(pk)) {
++idx;
reco::GenParticle & dummy_mother = cands->at(idx);
dummy_mother = reco::GenParticle(0, pk.p4(), pk.vertex(), 111, 2, true);
for (unsigned m = 0; m < pk.numberOfMothers(); ++m) {
new_cand.addMother(reco::GenParticleRef(ref, idx));
// Since the packed candidates drop the vertex position we'll take this from the mother
if (m == 0) {
dummy_mother.setP4(pk.mother(m)->p4());
dummy_mother.setVertex(pk.mother(m)->vertex());
new_cand.setVertex(pk.mother(m)->vertex());
}
// Should then add the GenParticle as a daughter of its dummy mother
dummy_mother.addDaughter(reco::GenParticleRef(ref, idx-1));
}
}
// Connect to mother from pruned particles
reco::GenParticle & daughter = cands->at(idx);
for (unsigned m = 0; m < pk.numberOfMothers(); ++m) {
daughter.addMother(reco::GenParticleRef(ref, pruned_idx_map.at(pk.mother(m))));
// Since the packed candidates drop the vertex position we'll take this from the mother
if (m == 0) {
daughter.setVertex(pk.mother(m)->vertex());
}
// Should then add this GenParticle as a daughter of its mother
cands->at(pruned_idx_map.at(pk.mother(m))).addDaughter(reco::GenParticleRef(ref, idx));
}
++idx;
}
event.put(std::move(cands));
}
bool MergedGenParticleProducer::isPhotonFromPrunedHadron(const pat::PackedGenParticle& pk) const
{
HepPDT::ParticleID motherid(pk.mother(0)->pdgId());
return
( pk.pdgId() == 22 // We care about photons for lepton dressing here
and pk.statusFlags().isDirectHadronDecayProduct() // Gen status flag seems correct
// Catch cases where miniaod mother is not compatible with the status flag
and not (motherid.isHadron() and pk.mother(0)->status() == 2)
);
}
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(MergedGenParticleProducer);