CMS 3D CMS Logo
GeneratorAnalyzer.py
Go to the documentation of this file.
1 from __future__ import print_function
2 from builtins import range
3 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
4 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
5 from PhysicsTools.Heppy.physicsutils.genutils import isNotFromHadronicShower, realGenMothers, realGenDaughters
6 
7 def interestingPdgId(id,includeLeptons=False):
8  id = abs(id)
9  return id in [6,7,8,17,18] or (includeLeptons and 11 <= id and id < 16) or (22 <= id and id < 40) or id > 1000000
10 
11 class GeneratorAnalyzer( Analyzer ):
12  """Save the hard-scattering final state of the event: top quarks, gauge & higgs bosons and BSM
13  particles, plus their immediate decay products, and their siblings (in order to get the jets
14  from matched X+jets generation.
15  Incoming partons are not included, by design.
16 
17  In the default configuration, leptons, light quarks and gluons are saved before FSR (a la status 3).
18  Everything else is saved after any radiation, i.e. immediately before the decay.
19 
20  Particles are saved in a list event.generatorSummary, with the index of their mothers ('motherIndex')
21  if the mother is also in the list, and with the pdgId of the mother ('motherId') and grand-mother
22  ('grandmotherId'). Particles also carry their index in the miniAOD genparticles collection ('rawIndex')
23  In addition, a 'sourceId' is set to the pdgId of the heaviest ancestor (or of the particle itself)
24  i.e. in top -> W -> lepton: the lepton sourceId will be 6
25  in tt+W with W -> lepton, the sourceId of the lepton will be 24.
26  sourceId will be 99 for paricles from hard scattering whose mother is light
27 
28  If requested, the full list of genParticles is also produced in event.genParticles (with indices
29  aligned to the miniAOD one). For particles that are in the generatorSummary, the same object is used.
30  An extra index 'genSummaryIndex' will be added to all particles, with the index in the generatorSummary
31  or -1 if the particle is not in the generatorSummary.
32 
33  Also, if requested it creates the splitted collections:
34  event.genHiggsBosons = []
35  event.genVBosons = []
36  event.gennus = [] # prompt neutrinos
37  event.gennusFromTop = [] # Neutrinos from t->W decay
38  event.genleps = [] # leptons from direct decays
39  event.gentauleps = [] # leptons from prompt taus
40  event.gentaus = [] # hadronically-decaying taus (if allGenTaus is False) or all taus (if allGenTaus is True)
41  event.gentopquarks = []
42  event.genbquarks = [] # b quarks from hard event (e.g. from top decays)
43  event.genwzquarks = [] # quarks from W,Z decays
44  event.genbquarksFromTop = []
45  event.genbquarksFromH = []
46  event.genlepsFromTop = [] #mu/ele that have a t->W chain as ancestor, also contained in event.genleps
47  event.genwzquarks and event.genbquarks, might have overlaps
48  event.genbquarksFromTop and event.genbquarksFromH are all contained in event.genbquarks
49 
50  """
51 
52  def __init__(self, cfg_ana, cfg_comp, looperName ):
53  super(GeneratorAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
54  self.stableBSMParticleIds = set(cfg_ana.stableBSMParticleIds) # neutralinos and such
55  self.savePreFSRParticleIds = set(cfg_ana.savePreFSRParticleIds)
56  self.makeAllGenParticles = cfg_ana.makeAllGenParticles
57  self.makeSplittedGenLists = cfg_ana.makeSplittedGenLists
58  self.allGenTaus = cfg_ana.allGenTaus if self.makeSplittedGenLists else False
59 
60  def declareHandles(self):
61  super(GeneratorAnalyzer, self).declareHandles()
62  self.mchandles['genParticles'] = AutoHandle( 'prunedGenParticles', 'std::vector<reco::GenParticle>' )
63 
64  def beginLoop(self,setup):
65  super(GeneratorAnalyzer,self).beginLoop(setup)
66 
67  def makeMCInfo(self, event):
68  verbose = getattr(self.cfg_ana, 'verbose', False)
69  rawGenParticles = self.mchandles['genParticles'].product()
70  good = []; keymap = {};
71  allGenParticles = []
72  for rawIndex,p in enumerate(rawGenParticles):
73  if self.makeAllGenParticles: allGenParticles.append(p)
74  id = abs(p.pdgId())
75  status = p.status()
76  # particles must be status > 2, except for prompt leptons, photons, neutralinos
77  if status <= 2:
78  if ((id not in self.stableBSMParticleIds) and
79  (id not in [11,12,13,14,15,16,22] or not isNotFromHadronicShower(p))):
80  continue
81  # a particle must not be decaying into itself
82  #print " test %6d : %+8d %3d : %8.2f %+5.2f %+5.2f : %d %d : %+8d {%6d}: %s" % ( rawIndex,
83  # p.pdgId(), p.status(), p.pt(), p.eta(), p.phi(), p.numberOfMothers(), p.numberOfDaughters(),
84  # p.motherRef(0).pdgId() if p.numberOfMothers() > 0 else -999, p.motherRef(0).key() if p.numberOfMothers() > 0 else -999,
85  # " ".join("%d[%d]" % (p.daughter(i).pdgId(), p.daughter(i).status()) for i in xrange(p.numberOfDaughters())))
86  if id in self.savePreFSRParticleIds:
87  # for light objects, we want them pre-radiation
88  if any((p.mother(j).pdgId() == p.pdgId()) for j in range(p.numberOfMothers())):
89  #print " fail auto-decay"
90  continue
91  else:
92  # everything else, we want it after radiation, i.e. just before decay
93  if any((p.daughter(j).pdgId() == p.pdgId() and p.daughter(j).status() > 2) for j in range(p.numberOfDaughters())):
94  #print " fail auto-decay"
95  continue
96  # FIXME find a better criterion to discard there
97  if status == 71:
98  #drop QCD radiation with unclear parentage
99  continue
100  # is it an interesting particle?
101  ok = False
102  if interestingPdgId(id):
103  #print " pass pdgId"
104  ok = True
105 
109  if not ok:
110  for mom in realGenMothers(p):
111  if interestingPdgId(mom.pdgId()) or (getattr(mom,'rawIndex',-1) in keymap):
112  #print " interesting mom"
113  # exclude extra x from p -> p + x
114  if not any(mom.daughter(j2).pdgId() == mom.pdgId() for j2 in range(mom.numberOfDaughters())):
115  #print " pass no-self-decay"
116  ok = True
117  # Account for generator feature with Higgs decaying to itself with same four-vector but no daughters
118  elif mom.pdgId() == 25 and any(mom.daughter(j2).pdgId() == 25 and mom.daughter(j2).numberOfDaughters()==0 for j2 in range(mom.numberOfDaughters())):
119  ok = True
120  if abs(mom.pdgId()) == 15:
121  # if we're a tau daughter we're status 2
122  # if we passed all the previous steps, then we're a prompt lepton
123  # so we'd like to be included
124  ok = True
125  if not ok and p.pt() > 10 and id in [1,2,3,4,5,21,22] and any(interestingPdgId(d.pdgId()) for d in realGenDaughters(mom)):
126  # interesting for being a parton brother of an interesting particle (to get the extra jets in ME+PS)
127  ok = True
128  if not ok and id in [11, 13, 15] and abs(mom.pdgId()) in [1,2,3,4,5,21]:
129  # Lepton e.g. in off-shell DY with the mother being one of the incoming partons
130  ok = True
131  if ok:
132  gp = p
133  gp.rawIndex = rawIndex # remember its index, so that we can set the mother index later
134  keymap[rawIndex] = len(good)
135  good.append(gp)
136  # connect mother links
137  for igp,gp in enumerate(good):
138  gp.motherIndex = -1
139  gp.sourceId = 99
140  gp.genSummaryIndex = igp
141  ancestor = None if gp.numberOfMothers() == 0 else gp.motherRef(0)
142  while ancestor != None and ancestor.isNonnull():
143  if ancestor.key() in keymap:
144  gp.motherIndex = keymap[ancestor.key()]
145  if ancestor.pdgId() != good[gp.motherIndex].pdgId():
146  print("Error keying %d: motherIndex %d, ancestor.pdgId %d, good[gp.motherIndex].pdgId() %d " % (igp, gp.motherIndex, ancestor.pdgId(), good[gp.motherIndex].pdgId()))
147  break
148  ancestor = None if ancestor.numberOfMothers() == 0 else ancestor.motherRef(0)
149  if abs(gp.pdgId()) not in [1,2,3,4,5,11,12,13,14,15,16,21]:
150  gp.sourceId = gp.pdgId()
151  if gp.motherIndex != -1:
152  ancestor = good[gp.motherIndex]
153  if ancestor.sourceId != 99 and (ancestor.mass() > gp.mass() or gp.sourceId == 99):
154  gp.sourceId = ancestor.sourceId
155  event.generatorSummary = good
156  # add the ID of the mother to be able to recreate last decay chains
157  for ip,p in enumerate(good):
158  moms = realGenMothers(p)
159  if len(moms)==0:
160  p.motherId = 0
161  p.grandmotherId = 0
162  elif len(moms)==1:
163  p.motherId = moms[0].pdgId()
164  gmoms = realGenMothers(moms[0])
165  p.grandmotherId = (gmoms[0].pdgId() if len(gmoms)==1 else (0 if len(gmoms)==0 else -9999))
166  else:
167  #print " unclear what mothers to give to this particle, among "," ".join("%d[%d]" % (m.pdgId(),m.status()) for m in moms)
168  p.motherId = -9999
169  p.grandmotherId = -9999
170  if verbose:
171  print("%3d {%6d}: %+8d %3d : %8.2f %+5.2f %+5.2f : %d %2d : %+8d {%3d}: %s" % ( ip,p.rawIndex,
172  p.pdgId(), p.status(), p.pt(), p.eta(), p.phi(), len(moms), p.numberOfDaughters(),
173  p.motherId, p.motherIndex,
174  " ".join("%d[%d]" % (p.daughter(i).pdgId(), p.daughter(i).status()) for i in range(p.numberOfDaughters()))))
175  if verbose:
176  print("\n\n")
177 
178  if self.makeAllGenParticles:
179  event.genParticles = allGenParticles
180 
181  if self.makeSplittedGenLists:
182  event.genHiggsBosons = []
183  event.genVBosons = []
184  event.gennus = []
185  event.gennusFromTop = []
186  event.genleps = []
187  event.gentauleps = []
188  event.gentaus = []
189  event.gentopquarks = []
190  event.genbquarks = []
191  event.genwzquarks = []
192  event.genbquarksFromTop = []
193  event.genbquarksFromH = []
194  event.genlepsFromTop = []
195  for p in event.generatorSummary:
196  id = abs(p.pdgId())
197  if id == 25:
198  event.genHiggsBosons.append(p)
199  elif id in {23,24}:
200  event.genVBosons.append(p)
201  elif id in {12,14,16}:
202  event.gennus.append(p)
203 
204  momids = [(m, abs(m.pdgId())) for m in realGenMothers(p)]
205 
206  #have a look at the lepton mothers
207  for mom, momid in momids:
208  #lepton from W
209  if momid == 24:
210  wmomids = [abs(m.pdgId()) for m in realGenMothers(mom)]
211  #W from t
212  if 6 in wmomids:
213  #save mu,e from t->W->mu/e
214  event.gennusFromTop.append(p)
215 
216  elif id in {11,13}:
217  #taus to separate vector
218  if abs(p.motherId) == 15:
219  event.gentauleps.append(p)
220  #all muons and electrons
221  else:
222  event.genleps.append(p)
223  momids = [(m, abs(m.pdgId())) for m in realGenMothers(p)]
224 
225  #have a look at the lepton mothers
226  for mom, momid in momids:
227  #lepton from W
228  if momid == 24:
229  wmomids = [abs(m.pdgId()) for m in realGenMothers(mom)]
230  #W from t
231  if 6 in wmomids:
232  #save mu,e from t->W->mu/e
233  event.genlepsFromTop.append(p)
234  elif id == 15:
235  if self.allGenTaus or not any([abs(d.pdgId()) in {11,13} for d in realGenDaughters(p)]):
236  event.gentaus.append(p)
237  elif id == 6:
238  event.gentopquarks.append(p)
239  elif id == 5:
240  event.genbquarks.append(p)
241  momids = [abs(m.pdgId()) for m in realGenMothers(p)]
242  if 6 in momids: event.genbquarksFromTop.append(p)
243  if 25 in momids: event.genbquarksFromH.append(p)
244  if id <= 5 and any([abs(m.pdgId()) in {23,24} for m in realGenMothers(p)]):
245  event.genwzquarks.append(p)
246 
247  def process(self, event):
248  self.readCollections( event.input )
249 
250  # if not MC, nothing to do
251  if not self.cfg_comp.isMC:
252  return True
253  # do MC level analysis
254  self.makeMCInfo(event)
255  return True
256 
257 import PhysicsTools.HeppyCore.framework.config as cfg
258 setattr(GeneratorAnalyzer,"defaultConfig",
259  cfg.Analyzer(GeneratorAnalyzer,
260  # BSM particles that can appear with status <= 2 and should be kept
261  stableBSMParticleIds = [ 1000022 ],
262  # Particles of which we want to save the pre-FSR momentum (a la status 3).
263  # Note that for quarks and gluons the post-FSR doesn't make sense,
264  # so those should always be in the list
265  savePreFSRParticleIds = [ 1,2,3,4,5, 11,12,13,14,15,16, 21 ],
266  # Make also the list of all genParticles, for other analyzers to handle
267  makeAllGenParticles = True,
268  # Make also the splitted lists
269  makeSplittedGenLists = True,
270  allGenTaus = False,
271  # Print out debug information
272  verbose = False,
273  )
274 )
FastTimerService_cff.range
range
Definition: FastTimerService_cff.py:34
gen.GeneratorAnalyzer.GeneratorAnalyzer.makeSplittedGenLists
makeSplittedGenLists
Definition: GeneratorAnalyzer.py:57
gen.GeneratorAnalyzer.interestingPdgId
def interestingPdgId(id, includeLeptons=False)
Definition: GeneratorAnalyzer.py:7
gen.GeneratorAnalyzer.GeneratorAnalyzer.process
def process(self, event)
Definition: GeneratorAnalyzer.py:247
mps_update.status
status
Definition: mps_update.py:69
join
static std::string join(char **cmd)
Definition: RemoteFile.cc:17
gen.GeneratorAnalyzer.GeneratorAnalyzer.makeAllGenParticles
makeAllGenParticles
Definition: GeneratorAnalyzer.py:56
gen.GeneratorAnalyzer.GeneratorAnalyzer.declareHandles
def declareHandles(self)
Definition: GeneratorAnalyzer.py:60
gen.GeneratorAnalyzer.GeneratorAnalyzer
Definition: GeneratorAnalyzer.py:11
genutils
any
bool any(const std::vector< T > &v, const T &what)
Definition: ECalSD.cc:38
gen.GeneratorAnalyzer.GeneratorAnalyzer.__init__
def __init__(self, cfg_ana, cfg_comp, looperName)
Definition: GeneratorAnalyzer.py:52
genutils.realGenMothers
def realGenMothers(gp)
Definition: genutils.py:100
gen.GeneratorAnalyzer.GeneratorAnalyzer.beginLoop
def beginLoop(self, setup)
Definition: GeneratorAnalyzer.py:64
remoteMonitoring_LED_IterMethod_cfg.Analyzer
Analyzer
Definition: remoteMonitoring_LED_IterMethod_cfg.py:32
print
void print(TMatrixD &m, const char *label=nullptr, bool mathematicaFormat=false)
Definition: Utilities.cc:46
EgammaValidation_cff.pdgId
pdgId
Definition: EgammaValidation_cff.py:118
gen.GeneratorAnalyzer.GeneratorAnalyzer.makeMCInfo
def makeMCInfo(self, event)
Definition: GeneratorAnalyzer.py:67
gen.GeneratorAnalyzer.GeneratorAnalyzer.allGenTaus
allGenTaus
Definition: GeneratorAnalyzer.py:58
genutils.isNotFromHadronicShower
def isNotFromHadronicShower(l)
Definition: genutils.py:66
gen.GeneratorAnalyzer.GeneratorAnalyzer.stableBSMParticleIds
stableBSMParticleIds
Definition: GeneratorAnalyzer.py:54
gen.GeneratorAnalyzer.GeneratorAnalyzer.savePreFSRParticleIds
savePreFSRParticleIds
Definition: GeneratorAnalyzer.py:55
genutils.realGenDaughters
def realGenDaughters(gp, excludeRadiation=True)
Definition: genutils.py:81
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Generated for CMSSW Reference Manual by doxygen 1.8.16