analyzers/objects/VertexAnalyzer.py

0001 from __future__ import print_function
0002 import itertools
0003
0004 from PhysicsTools.Heppy.analyzers.core.VertexHistograms import VertexHistograms
0005 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
0006 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
0007 from PhysicsTools.HeppyCore.statistics.average import Average
0008 from PhysicsTools.Heppy.physicsutils.PileUpSummaryInfo import PileUpSummaryInfo
0009 import PhysicsTools.HeppyCore.framework.config as cfg
0010
0011 class VertexAnalyzer( Analyzer ):
0012     """selects a list of good primary vertices,
0013     and optionally add a pile-up weight to MC events.
0014
0015     The list of good primary vertices is put in event.goodVertices.
0016     if no good vertex is found, the process function returns False.
0017
0018     The weight is put in event.vertexWeight, and is multiplied to
0019     the global event weight, event.eventWeight.
0020
0021     Example:
0022
0023     vertexAna = cfg.Analyzer(
0024       'VertexAnalyzer',
0025       goodVertices = 'goodPVFilter',
0026       vertexWeight = 'vertexWeightFall112011AB',
0027       # uncomment the following line if you want a vertex weight = 1 (no weighting)
0028       # fixedWeight = 1,
0029       verbose = False
0030       )
0031
0032     If fixedWeight is set to None, the vertex weight is read from the EDM collection with module name
0033     'vertexWeightFall112011AB'.
0034     Otherwise, the weight is set to fixedWeight.
0035
0036     The vertex weight collection was at some point produced in the PAT+CMG step,
0037     and could directly be accessed from the PAT or CMG tuple.
0038     In the most recent versions of the PAT+CMG tuple, this collection is not present anymore,
0039     and an additional full framework process must be ran to produce this collection,
0040     so that this analyzer can read it. An example cfg to do that can be found here:
0041     http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/CMG/CMGTools/H2TauTau/prod/vertexWeight2011_cfg.py?view=markup
0042
0043
0044     """
0045
0046     def __init__(self, cfg_ana, cfg_comp, looperName):
0047         super(VertexAnalyzer, self).__init__(cfg_ana, cfg_comp, looperName)
0048
0049         self.doHists=True
0050         if (hasattr(self.cfg_ana,'makeHists')) and (not self.cfg_ana.makeHists):
0051             self.doHists=False
0052         if self.doHists:
0053             self.pileup = VertexHistograms('/'.join([self.dirName,
0054                                                      'pileup.root']))
0055
0056         self.allVertices = self.cfg_ana.allVertices if (hasattr(self.cfg_ana,'allVertices')) else "_AUTO_"
0057
0058     def declareHandles(self):
0059         super(VertexAnalyzer, self).declareHandles()
0060         if self.allVertices == '_AUTO_':
0061           self.handles['vertices'] =  AutoHandle( "offlineSlimmedPrimaryVertices", 'std::vector<reco::Vertex>', fallbackLabel="offlinePrimaryVertices" )
0062         else:
0063           self.handles['vertices'] =  AutoHandle( self.allVertices, 'std::vector<reco::Vertex>' )
0064         self.fixedWeight = None
0065         if self.cfg_comp.isMC:
0066             if hasattr( self.cfg_ana, 'fixedWeight'):
0067                 self.fixedWeight = self.cfg_ana.fixedWeight
0068             else:
0069                 self.mchandles['vertexWeight'] = AutoHandle( self.cfg_ana.vertexWeight,
0070                                                              'double' )
0071
0072         self.mchandles['pusi'] =  AutoHandle(
0073             'slimmedAddPileupInfo',
0074             'std::vector<PileupSummaryInfo>',
0075             fallbackLabel='addPileupInfo',
0076             )
0077
0078         self.handles['rho'] =  AutoHandle(
0079             ('fixedGridRhoFastjetAll',''),
0080             'double'
0081             )
0082         self.handles['rhoCN'] =  AutoHandle(
0083             ('fixedGridRhoFastjetCentralNeutral',''),
0084             'double'
0085             )
0086         self.handles['sigma'] =  AutoHandle(
0087             ('fixedGridSigmaFastjetAll',''),
0088             'double',
0089             mayFail=True
0090             )
0091
0092     def beginLoop(self, setup):
0093         super(VertexAnalyzer,self).beginLoop(setup)
0094         self.averages.add('vertexWeight', Average('vertexWeight') )
0095         self.counters.addCounter('GoodVertex')
0096         self.count = self.counters.counter('GoodVertex')
0097         self.count.register('All Events')
0098         self.count.register('Events With Good Vertex')
0099
0100
0101     def process(self,  event):
0102         self.readCollections(event.input )
0103         event.rho = self.handles['rho'].product()[0]
0104         event.rhoCN = self.handles['rhoCN'].product()[0]
0105         event.sigma = self.handles['sigma'].product()[0] if self.handles['sigma'].isValid() else -999
0106         event.vertices = self.handles['vertices'].product()
0107         event.goodVertices = list(filter(self.testGoodVertex,event.vertices))
0108
0109
0110         self.count.inc('All Events')
0111
0112
0113         event.vertexWeight = 1
0114         if self.cfg_comp.isMC:
0115             event.pileUpInfo = map( PileUpSummaryInfo,
0116                                     self.mchandles['pusi'].product() )
0117             if self.fixedWeight is None:
0118                 event.vertexWeight = self.mchandles['vertexWeight'].product()[0]
0119             else:
0120                 event.vertexWeight = self.fixedWeight
0121         event.eventWeight *= event.vertexWeight
0122
0123         self.averages['vertexWeight'].add( event.vertexWeight )
0124         if self.verbose:
0125             print('VertexAnalyzer: #vert = ', len(event.vertices), \
0126                   ', weight = ', event.vertexWeight)
0127
0128         # Check if events needs to be skipped if no good vertex is found (useful for generator level studies)
0129         keepFailingEvents = False
0130         if hasattr( self.cfg_ana, 'keepFailingEvents'):
0131             keepFailingEvents = self.cfg_ana.keepFailingEvents
0132         if len(event.goodVertices)==0:
0133             event.passedVertexAnalyzer=False
0134             if not keepFailingEvents:
0135                 return False
0136         else:
0137             event.passedVertexAnalyzer=True
0138
0139         if self.doHists:
0140             self.pileup.hist.Fill( len(event.goodVertices) )
0141 #A.R. mindist is one of the slowest functions, default commented
0142 #           self.pileup.mindist.Fill( self.mindist(event.goodVertices) )
0143
0144         self.count.inc('Events With Good Vertex')
0145         return True
0146
0147
0148     def testGoodVertex(self,vertex):
0149         if vertex.isFake():
0150             return False
0151         if vertex.ndof()<=4:
0152             return False
0153         if abs(vertex.z())>24:
0154             return False
0155         if vertex.position().Rho()>2:
0156             return False
0157
0158         return True
0159
0160     def mindist(self, vertices):
0161         mindist = 999999
0162         for comb in itertools.combinations(vertices, 2):
0163             dist = abs(comb[0].z() - comb[1].z())
0164             if dist<mindist:
0165                 mindist = dist
0166         return mindist
0167
0168     def write(self, setup):
0169         super(VertexAnalyzer, self).write(setup)
0170         if self.doHists:
0171             self.pileup.write()
0172
0173 setattr(VertexAnalyzer,"defaultConfig",cfg.Analyzer(
0174     class_object=VertexAnalyzer,
0175     vertexWeight = None,
0176     fixedWeight = 1,
0177     verbose = False
0178    )
0179 )