analyzers/objects/VertexAnalyzer.py

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