diff --git a/CMGTools/RootTools/python/analyzers/skimAnalyzerCount.py b/CMGTools/RootTools/python/analyzers/skimAnalyzerCount.py index 75d8bc560e161..edf7cd36189a0 100644 --- a/CMGTools/RootTools/python/analyzers/skimAnalyzerCount.py +++ b/CMGTools/RootTools/python/analyzers/skimAnalyzerCount.py @@ -14,7 +14,7 @@ class skimAnalyzerCount( Analyzer ): def __init__(self, cfg_ana, cfg_comp, looperName): super(skimAnalyzerCount, self).__init__(cfg_ana, cfg_comp, looperName) - self.useLumiBlocks = self.cfg_ana.useLumiBlocks if (hasattr(self.cfg_ana,'useLumiBlocks')) else True + self.useLumiBlocks = self.cfg_ana.useLumiBlocks if (hasattr(self.cfg_ana,'useLumiBlocks')) else False def declareHandles(self): super(skimAnalyzerCount, self).declareHandles() @@ -27,7 +27,8 @@ def beginLoop(self): self.counters.addCounter('SkimReport') self.count = self.counters.counter('SkimReport') self.count.register('All Events') - self.count.register('Sum Weights') + if self.cfg_comp.isMC: + self.count.register('Sum Weights') if not self.useLumiBlocks: print 'Will actually count events instead of accessing lumi blocks' @@ -47,7 +48,8 @@ def beginLoop(self): if self.useLumiBlocks: self.count.inc('All Events',totalEvents) - self.count.inc('Sum Weights',totalEvents) + if self.cfg_comp.isMC: + self.count.inc('Sum Weights',totalEvents) print 'Done -> proceeding with the analysis' else: print 'Failed -> will have to actually count events (this can happen if the input dataset is not a CMG one)' @@ -58,5 +60,6 @@ def process(self, iEvent, event): if not self.useLumiBlocks: self.readCollections( iEvent ) self.count.inc('All Events') - self.count.inc('Sum Weights', self.mchandles['GenInfo'].product().weight()) + if self.cfg_comp.isMC: + self.count.inc('Sum Weights', self.mchandles['GenInfo'].product().weight()) return True diff --git a/CMGTools/RootTools/python/fwlite/Looper.py b/CMGTools/RootTools/python/fwlite/Looper.py index 07bdca202aadc..9339d078fc06e 100644 --- a/CMGTools/RootTools/python/fwlite/Looper.py +++ b/CMGTools/RootTools/python/fwlite/Looper.py @@ -4,6 +4,7 @@ # import glob import logging import pprint +import time from DataFormats.FWLite import Events, Handle from CMGTools.RootTools.fwlite.Event import Event from CMGTools.RootTools.fwlite.PythonPath import pythonpath @@ -35,7 +36,7 @@ def __init__( self, name, cfg_comp, sequence, nEvents=None, firstEvent=0, nPrint self.classes = {} #TODO: should be a diclist? self.analyzers = map( self._buildAnalyzer, sequence ) - self.nEvents = nEvents + self.nEvents = getattr(cfg_comp, 'maxEvents', nEvents) self.firstEvent = firstEvent self.nPrint = int(nPrint) # initialize FWLite chain on input file: @@ -118,7 +119,11 @@ def loop(self): # if iEv == nEvents: # break if iEv%100 ==0: - print 'event', iEv + if iEv == 100: + print 'event', iEv + self.start_time = time.time() + elif iEv > 100: + print 'event %d (%.1f ev/s)' % (iEv, (iEv-100)/float(time.time() - self.start_time)) self.process( iEv ) if iEv= 2 loose leptons, no pt cuts or id ## +########################################################## + +import CMGTools.RootTools.fwlite.Config as cfg +from CMGTools.RootTools.fwlite.Config import printComps +from CMGTools.RootTools.RootTools import * + +#Load all analyzers +from CMGTools.TTHAnalysis.analyzers.susyCore_modules_cff import * + +# --- LEPTON SKIMMING --- +ttHLepSkim.minLeptons = 2 +ttHLepSkim.maxLeptons = 999 + +ttHGenLevel = cfg.Analyzer( + 'ttHGenLevelOnlyStudy', + muon_pt_min = 5., + electron_pt_min = 7., +) + + +from CMGTools.TTHAnalysis.samples.samples_8TeV_v517 import triggers_mumu, triggers_ee, triggers_mue, triggers_1mu +# Tree Producer +treeProducer = cfg.Analyzer( + 'treeProducerSusyGenLevelOnly', + vectorTree = True, + saveTLorentzVectors = False, # can set to True to get also the TLorentzVectors, but trees will be bigger + PDFWeights = PDFWeights, + # triggerBits = {} # no HLT + ) + + +#-------- SAMPLES AND TRIGGERS ----------- +from CMGTools.TTHAnalysis.samples.samples_8TeV_v517 import * +Test = kreator.makePrivateMCComponent('Test', '/store/cmst3/user/gpetrucc/maiani', ['m100_g050_3mu.GEN.root'] ) +WZ3l_ascms = kreator.makePrivateMCComponent('WZ3l_ascms', '/store/cmst3/user/gpetrucc/maiani/tests', ['xs_wz_3l_ascms.GEN.root']) +WZ3mu_ascms = kreator.makePrivateMCComponent('WZ3mu_ascms', '/store/cmst3/user/gpetrucc/maiani/tests', ['xs_wz_3mu_ascms.GEN.root']) +WZ3mu_offshell = kreator.makePrivateMCComponent('WZ3mu_offshell', '/store/cmst3/user/gpetrucc/maiani/tests', ['xs_wz_3mu_offshell.GEN.root']) + +GEN_S3m_lo = kreator.makePrivateMCComponent('GEN_S3m_lo', '/store/cmst3/user/gpetrucc/lmutau/madtests/', [ 'lmutau_lo_test.GEN.root' ]) +GEN_S3m_lo_012j = kreator.makePrivateMCComponent('GEN_S3m_lo_012j', '/store/cmst3/user/gpetrucc/lmutau/madtests/', [ 'lmutau_lo_012j_test.GEN.root' ]) +GEN_S3m_lo_01j = kreator.makePrivateMCComponent('GEN_S3m_lo_01j', '/store/cmst3/user/gpetrucc/lmutau/madtests/', [ 'lmutau_lo_01j_test.GEN.root' ]) +GEN_S3m_nlo = kreator.makePrivateMCComponent('GEN_S3m_nlo', '/store/cmst3/user/gpetrucc/lmutau/madtests/', [ 'lmutau_nlo_test.GEN.root' ]) +GEN_S3m_nlo_01j = kreator.makePrivateMCComponent('GEN_S3m_nlo_01j', '/store/cmst3/user/gpetrucc/lmutau/madtests/', [ 'lmutau_nlo_01j_test.GEN.root', 'lmutau_nlo_01j_test.2.GEN.root' ]) +GEN_S3m_lo_direct = kreator.makePrivateMCComponent('GEN_S3m_lo_direct', '/store/cmst3/user/gpetrucc/lmutau/madtests/', [ 'lmutau_lo_direct_test.GEN.root' ]) + +GEN_Bonly_S3m = kreator.makePrivateMCComponent('GEN_Bonly_S3m', '/store/cmst3/user/gpetrucc/lmutau/', [ 'Bonly_S3m_m100_g050_8TeV.GEN.root' ]) +GEN_SBI_S3m_m100_g050 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m100_g050', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m100_g050_8TeV.GEN.root' ]) +GEN_SBI_S3m_m105_g052 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m105_g052', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m105_g052_8TeV.GEN.root' ]) +GEN_SBI_S3m_m110_g055 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m110_g055', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m110_g055_8TeV.GEN.root' ]) +GEN_SBI_S3m_m115_g057 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m115_g057', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m115_g057_8TeV.GEN.root' ]) +GEN_SBI_S3m_m120_g060 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m120_g060', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m120_g060_8TeV.GEN.root' ]) +GEN_SBI_S3m_m125_g062 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m125_g062', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m125_g062_8TeV.GEN.root' ]) +GEN_SBI_S3m_m70_g035 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m70_g035', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m70_g035_8TeV.GEN.root' ]) +GEN_SBI_S3m_m75_g037 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m75_g037', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m75_g037_8TeV.GEN.root' ]) +GEN_SBI_S3m_m80_g040 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m80_g040', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m80_g040_8TeV.GEN.root' ]) +GEN_SBI_S3m_m85_g042 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m85_g042', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m85_g042_8TeV.GEN.root' ]) +GEN_SBI_S3m_m90_g045 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m90_g045', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m90_g045_8TeV.GEN.root' ]) +GEN_SBI_S3m_m95_g047 = kreator.makePrivateMCComponent('GEN_SBI_S3m_m95_g047', '/store/cmst3/user/gpetrucc/lmutau/', [ 'SBI_S3m_m95_g047_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m100_g050 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m100_g050', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m100_g050_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m105_g052 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m105_g052', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m105_g052_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m110_g055 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m110_g055', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m110_g055_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m115_g057 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m115_g057', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m115_g057_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m120_g060 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m120_g060', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m120_g060_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m125_g062 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m125_g062', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m125_g062_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m70_g035 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m70_g035', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m70_g035_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m75_g037 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m75_g037', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m75_g037_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m80_g040 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m80_g040', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m80_g040_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m85_g042 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m85_g042', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m85_g042_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m90_g045 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m90_g045', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m90_g045_8TeV.GEN.root' ]) +GEN_Sonly_S3m_m95_g047 = kreator.makePrivateMCComponent('GEN_Sonly_S3m_m95_g047', '/store/cmst3/user/gpetrucc/lmutau/v1/', [ 'S3m_m95_g047_8TeV.GEN.root' ]) + + +### ====== SUSY: PRIVATE PRODUCTIONS ========== +GEN_T2tt_py8had = kreator.makePrivateMCComponent('GEN_T2tt_py8had', '/store/cmst3/user/gpetrucc/SUSY/TestProd/T2tt/', [ 'T2tt_onshell_pyt8had.root' ] + [ "T2tt_onshell_pyt8had.run_0%d_chunk_%d.root" % (i,j) for i in 2,3,4,5, for j in 0,1,2,3,4 ] ) +GEN_T2tt_py8dec_py8had = kreator.makePrivateMCComponent('GEN_T2tt_py8dec_py8had', '/store/cmst3/user/gpetrucc/SUSY/TestProd/T2tt/', [ 'T2tt_onshell_py8decay_pyt8had.root' ] + [ "T2tt_onshell_py8decay_pyt8had.run_0%d_chunk_%d.root" % (i,j) for i in 2,3,4,5, for j in 0,1,2,3,4 ]) +GEN_T2tt_mgdec_py8had = kreator.makePrivateMCComponent('GEN_T2tt_mgdec_py8had', '/store/cmst3/user/gpetrucc/SUSY/TestProd/T2tt/', [ 'T2tt_decayed_pyt8had.root' ] + [ "T2tt_decayed_pyt8had.run_0%d_chunk_%d.root" % (i,j) for i in 3,4,5,6, for j in 0,1,2,3,4 ]) + +GEN_T2tt_py8had_ch = kreator.makePrivateMCComponent('GEN_T2tt_py8had_ch', '/store/cmst3/user/gpetrucc/SUSY/TestProd/T2tt/', [ 'T2tt_onshell_pyt8had_chargino.root' ]+[ "T2tt_onshell_pyt8had_chargino.run_0%d_chunk_0%d.root" % (i,j) for i in 2,3,4,5, for j in 0,1,2,3,4 ] ) +GEN_T2tt_mgdec_py8had_ch = kreator.makePrivateMCComponent('GEN_T2tt_mgdec_py8had_ch', '/store/cmst3/user/gpetrucc/SUSY/TestProd/T2tt/', [ "T2tt_decayed_pyt8had_chargino.run_0%d_chunk_0%d.root" % (i,j) for i in (1,) for j in 0,1,2,3,4 ]) +GEN_T2tt_mgdec_py8had_both = kreator.makePrivateMCComponent('GEN_T2tt_mgdec_py8had_both', '/store/cmst3/user/gpetrucc/SUSY/TestProd/T2tt/', [ "T2tt_decayed_pyt8had_both.run_0%d_chunk_0%d.root" % (i,j) for i in 1,2 for j in 0,1,2,3,4 ]) + +GEN_T1tttt_mGo800_mStop300_mChi280_mg5 = kreator.makePrivateMCComponent('GEN_T1tttt_mGo800_mStop300_mChi280_mg5', '/store/cmst3/user/gpetrucc/SUSY/Prod/T1tttt_mGo800_mStop300_mChi280_mg5dec_pythia8/', [ "T1tttt_mGo800_mStop300_mChi280_mg5dec_pythia8.run_%02d_chunk_%02d.root" % (i,j) for i in (1,2) for j in xrange(10) ]) +#GEN_T1tttt_mGo1300_mStop300_mChi280_mg5 = kreator.makePrivateMCComponent('GEN_T1tttt_mGo1300_mStop300_mChi280_mg5', '/store/cmst3/user/gpetrucc/SUSY/Prod/T1tttt_mGo1300_mStop300_mChi280_mg5dec_pythia8/', [ "T1tttt_mGo1300_mStop300_mChi280_mg5dec_pythia8.run_%02d_chunk_%02d.root" % (i,j) for i in (1,2) for j in xrange(10) ]) +GEN_T1tttt_mGo800_mStop300_mCh285_mChi280_mg5 = kreator.makePrivateMCComponent('GEN_T1tttt_mGo800_mStop300_mCh285_mChi280_mg5', '/store/cmst3/user/gpetrucc/SUSY/Prod/T1tttt_mGo800_mStop300_mCh285_mChi280_mg5dec_pythia8/', [ "T1tttt_mGo800_mStop300_mCh285_mChi280_mg5dec_pythia8.run_%02d_chunk_%d.root" % (i,j) for i in (1,) for j in xrange(20) ]) +#GEN_T1tttt_mGo1300_mStop300_mCh285_mChi280_mg5 = kreator.makePrivateMCComponent('GEN_T1tttt_mGo1300_mStop300_mCh285_mChi280_mg5', '/store/cmst3/user/gpetrucc/SUSY/Prod/T1tttt_mGo1300_mStop300_mCh285_mChi280_mg5dec_pythia8/', [ "T1tttt_mGo1300_mStop300_mCh285_mChi280_mg5dec_pythia8.run_%02d_chunk_%02d.root" % (i,j) for i in (1,2) for j in xrange(10) ]) + +### ====== SUSY: PRIVATE DECAY+HADRONIZATION OF CENTRALLY PRODUCED LHE FILES ========== +eosGenFiles = [ x.strip() for x in open(os.environ["CMSSW_BASE"]+"/src/CMGTools/TTHAnalysis/python/samples/genLevel-susySMS-13TeV", "r") ] +print eosGenFiles +def mkGen(name,what): + return kreator.makePrivateMCComponent(name, '/'+what, [ f for f in eosGenFiles if ("/%s/"%what) in f ] ) +GEN_T1tttt_2J_mGo800_mStop300_mChi280_py8 = mkGen('GEN_T1tttt_2J_mGo800_mStop300_mChi280_py8', 'T1ttt_2J_mGo800_mStop300_mChi280_pythia8-4bodydec') +GEN_T1tttt_2J_mGo1300_mStop300_mChi280_py8 = mkGen('GEN_T1tttt_2J_mGo1300_mStop300_mChi280_py8', 'T1ttt_2J_mGo1300_mStop300_mChi280_pythia8-4bodydec') +GEN_T1tttt_2J_mGo800_mStop300_mCh285_mChi280_py8 = mkGen('GEN_T1tttt_2J_mGo800_mStop300_mCh285_mChi280_py8', 'T1ttt_2J_mGo800_mStop300_mCh285_mChi280_pythia8-23bodydec') +GEN_T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_py8 = mkGen('GEN_T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_py8', 'T1ttt_2J_mGo1300_mStop300_mCh285_mChi280_pythia8-23bodydec') +GEN_T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_py8_dilep = mkGen('GEN_T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_py8_dilep', 'T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_23bodydec_dilepfilter') + + +#selectedComponents = [ GEN_S3m_lo, GEN_S3m_lo_012j, GEN_S3m_lo_01j, GEN_S3m_nlo, GEN_S3m_nlo_01j, GEN_S3m_lo_direct ] +selectedComponents = [ GEN_T2tt_py8had, GEN_T2tt_py8dec_py8had, GEN_T2tt_mgdec_py8had ] +#selectedComponents = [ GEN_Bonly_S3m, GEN_SBI_S3m_m100_g050, GEN_SBI_S3m_m105_g052, GEN_SBI_S3m_m110_g055, GEN_SBI_S3m_m115_g057, GEN_SBI_S3m_m120_g060, GEN_SBI_S3m_m125_g062, GEN_SBI_S3m_m70_g035, GEN_SBI_S3m_m75_g037, GEN_SBI_S3m_m80_g040, GEN_SBI_S3m_m85_g042, GEN_SBI_S3m_m90_g045, GEN_SBI_S3m_m95_g047, GEN_Sonly_S3m_m100_g050, GEN_Sonly_S3m_m105_g052, GEN_Sonly_S3m_m110_g055, GEN_Sonly_S3m_m115_g057, GEN_Sonly_S3m_m120_g060, GEN_Sonly_S3m_m125_g062, GEN_Sonly_S3m_m70_g035, GEN_Sonly_S3m_m75_g037, GEN_Sonly_S3m_m80_g040, GEN_Sonly_S3m_m85_g042, GEN_Sonly_S3m_m90_g045, GEN_Sonly_S3m_m95_g047 ] +#selectedComponents = [ GEN_Bonly_S3m, GEN_SBI_S3m_m110_g055, GEN_Sonly_S3m_m110_g055 ] +selectedComponents = [ GEN_T1tttt_2J_mGo800_mStop300_mChi280_py8, GEN_T1tttt_2J_mGo1300_mStop300_mChi280_py8, GEN_T1tttt_2J_mGo800_mStop300_mCh285_mChi280_py8, GEN_T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_py8 ] +#selectedComponents = [ GEN_T2tt_mgdec_py8had_ch, GEN_T2tt_mgdec_py8had_both, GEN_T2tt_py8had_ch ] +#selectedComponents = [ GEN_T2tt_py8had, GEN_T2tt_mgdec_py8had ] +selectedComponents = [ GEN_T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_py8, GEN_T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_py8_dilep ] + +for c in selectedComponents: c.splitFactor = 100 +#-------- SEQUENCE + +sequence = cfg.Sequence([ + skimAnalyzer, + ttHGenLevel, + ttHLepSkim, + treeProducer, + ]) + + +#-------- HOW TO RUN +test = 0 +if test==1: + # test a single component, using a single thread. + comp = GEN_T2tt_mgdec_py8had_both + comp.files = comp.files[:1] + selectedComponents = [comp] + comp.splitFactor = 1 +elif test==2: + # test all components (1 thread per component). + for comp in selectedComponents: + comp.splitFactor = 1 + comp.files = comp.files[:1] + + + +config = cfg.Config( components = selectedComponents, + sequence = sequence ) + +printComps(config.components, True) diff --git a/CMGTools/TTHAnalysis/python/analyzers/ntupleTypes.py b/CMGTools/TTHAnalysis/python/analyzers/ntupleTypes.py index e14d69fe2be85..549a1384540b3 100644 --- a/CMGTools/TTHAnalysis/python/analyzers/ntupleTypes.py +++ b/CMGTools/TTHAnalysis/python/analyzers/ntupleTypes.py @@ -284,6 +284,7 @@ NTupleVariable("eip3d", lambda x : x.d3d.error(), help="Uncertainty on the 3D distance from the PV [cm]"), NTupleVariable("sip3d", lambda x : x.d3d.significance(), help="S_{ip3d} with respect to PV (absolute value)"), NTupleVariable("cosTheta", lambda x : x.cosTheta, help="Cosine of the angle between the 3D displacement and the momentum"), + NTupleVariable("mva", lambda x : x.mva, help="MVA discriminator"), NTupleVariable("jetPt", lambda x : x.jet.pt() if x.jet != None else 0, help="pT of associated jet"), NTupleVariable("jetBTag", lambda x : x.jet.btag('combinedSecondaryVertexBJetTags') if x.jet != None else -99, help="CSV b-tag of associated jet"), NTupleVariable("mcMatchNTracks", lambda x : x.mcMatchNTracks, int, mcOnly=True, help="Number of mc-matched tracks in SV"), diff --git a/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyGenLevelOnly.py b/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyGenLevelOnly.py new file mode 100644 index 0000000000000..c63e2f9a81ef8 --- /dev/null +++ b/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyGenLevelOnly.py @@ -0,0 +1,73 @@ +from CMGTools.TTHAnalysis.analyzers.ttHLepTreeProducerNew import * + +class treeProducerSusyGenLevelOnly( ttHLepTreeProducerNew ): + + #----------------------------------- + # CORE TREE PRODUCER FOR THE SUSY ANALYSES + # defines the core variables that will be present in the trees of all final states + #----------------------------------- + def __init__(self, cfg_ana, cfg_comp, looperName): + super(treeProducerSusyGenLevelOnly,self).__init__(cfg_ana, cfg_comp, looperName) + + ## Declare what we want to fill + self.globalVariables = [ + #NTupleVariable("rho", lambda ev: ev.rho, float, help="kt6PFJets rho"), + #NTupleVariable("nVert", lambda ev: len(ev.goodVertices), int, help="Number of good vertices"), + NTupleVariable("nJet25", lambda ev: len(ev.cleanJets), int, help="Number of jets with pt > 25"), + NTupleVariable("nJet40", lambda ev: sum([j.pt() > 40 for j in ev.cleanJets]), int, help="Number of jets with pt > 40, |eta|<2.4"), + NTupleVariable("nJet40a", lambda ev: sum([j.pt() > 40 for j in ev.cleanJetsAll]), int, help="Number of jets with pt > 40, |eta|<4.7"), + NTupleVariable("nBJetLoose25", lambda ev: sum([j.btagWP("CSVL") for j in ev.cleanJets if j.pt() > 25]), int, help="Number of jets with pt > 25 passing CSV loose"), + NTupleVariable("nBJetMedium25", lambda ev: sum([j.btagWP("CSVM") for j in ev.cleanJets if j.pt() > 25]), int, help="Number of jets with pt > 25 passing CSV medium"), + NTupleVariable("nBJetLoose40", lambda ev: sum([j.btagWP("CSVL") for j in ev.cleanJets if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CSV loose"), + NTupleVariable("nBJetMedium40", lambda ev: sum([j.btagWP("CSVM") for j in ev.cleanJets if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CSV medium"), + ##-------------------------------------------------- + NTupleVariable("nLepGood20", lambda ev: sum([l.pt() > 20 for l in ev.selectedLeptons]), int, help="Number of leptons with pt > 20"), + NTupleVariable("nLepGood15", lambda ev: sum([l.pt() > 15 for l in ev.selectedLeptons]), int, help="Number of leptons with pt > 15"), + NTupleVariable("nLepGood10", lambda ev: sum([l.pt() > 10 for l in ev.selectedLeptons]), int, help="Number of leptons with pt > 10"), + ##-------- custom jets ------------------------------------------ + NTupleVariable("htJet25", lambda ev : ev.htJet25, help="H_{T} computed from leptons and jets (with |eta|<2.4, pt > 25 GeV)"), + NTupleVariable("mhtJet25", lambda ev : ev.mhtJet25, help="H_{T}^{miss} computed from leptons and jets (with |eta|<2.4, pt > 25 GeV)"), + NTupleVariable("htJet40j", lambda ev : ev.htJet40j, help="H_{T} computed from only jets (with |eta|<2.4, pt > 40 GeV)"), + #NTupleVariable("htJet40ja", lambda ev : ev.htJet40j, help="H_{T} computed from only jets (with |eta|<4.7, pt > 40 GeV)"), + NTupleVariable("htJet40", lambda ev : ev.htJet40, help="H_{T} computed from leptons and jets (with |eta|<2.4, pt > 40 GeV)"), + #NTupleVariable("htJet40a", lambda ev : ev.htJet40a, help="H_{T} computed from leptons and jets (with |eta|<4.7, pt > 40 GeV)"), + NTupleVariable("mhtJet40", lambda ev : ev.mhtJet40, help="H_{T}^{miss} computed from leptons and jets (with |eta|<2.4, pt > 40 GeV)"), + #NTupleVariable("mhtJet40a", lambda ev : ev.mhtJet40a, help="H_{T}^{miss} computed from leptons and jets (with |eta|<4.7, pt > 40 GeV)"), + ##-------------------------------------------------- + NTupleVariable("mZ1", lambda ev : ev.bestZ1[0], help="Best m(ll) SF/OS"), + NTupleVariable("mZ1SFSS", lambda ev : ev.bestZ1sfss[0], help="Best m(ll) SF/SS"), + NTupleVariable("minMllSFOS", lambda ev: ev.minMllSFOS, help="min m(ll), SF/OS"), + NTupleVariable("maxMllSFOS", lambda ev: ev.maxMllSFOS, help="max m(ll), SF/OS"), + NTupleVariable("minMllAFOS", lambda ev: ev.minMllAFOS, help="min m(ll), AF/OS"), + NTupleVariable("maxMllAFOS", lambda ev: ev.maxMllAFOS, help="max m(ll), AF/OS"), + NTupleVariable("minMllAFSS", lambda ev: ev.minMllAFSS, help="min m(ll), AF/SS"), + NTupleVariable("maxMllAFSS", lambda ev: ev.maxMllAFSS, help="max m(ll), AF/SS"), + NTupleVariable("minMllAFAS", lambda ev: ev.minMllAFAS, help="min m(ll), AF/AS"), + NTupleVariable("maxMllAFAS", lambda ev: ev.maxMllAFAS, help="max m(ll), AF/AS"), + NTupleVariable("m2l", lambda ev: ev.m2l, help="m(ll)"), + ##-------------------------------------------------- + NTupleVariable("mZ2", lambda ev : ev.bestZ2[3], help="m(ll) of second SF/OS pair, for ZZ reco."), + NTupleVariable("m3l", lambda ev: ev.m3l, help="m(3l)"), + NTupleVariable("m4l", lambda ev: ev.m4l, help="m(4l)"), + NTupleVariable("pt2l", lambda ev: ev.pt2l, help="p_{T}(ll)"), + NTupleVariable("pt3l", lambda ev: ev.pt3l, help="p_{T}(3l)"), + NTupleVariable("pt4l", lambda ev: ev.pt4l, help="p_{T}(4l)"), + NTupleVariable("ht3l", lambda ev: ev.ht3l, help="H_{T}(3l)"), + NTupleVariable("ht4l", lambda ev: ev.ht4l, help="H_{T}(4l)"), + NTupleVariable("q3l", lambda ev: ev.q3l, int, help="q(3l)"), + NTupleVariable("q4l", lambda ev: ev.q4l, int, help="q(4l)"), + ] + + self.globalObjects = { + "met" : NTupleObject("met", fourVectorType, help="PF E_{T}^{miss}, after type 1 corrections"), + } + self.collections = { + "selectedLeptons" : NTupleCollection("LepGood", particleType, 8, help="Leptons after the preselection"), + "cleanJets" : NTupleCollection("Jet", fourVectorType, 8, help="Cental jets after full selection and cleaning, sorted by pt"), + "cleanJetsFwd" : NTupleCollection("JetFwd", fourVectorType, 4, help="Forward jets after full selection and cleaning, sorted by pt"), + } + + ## Book the variables, but only if we're called explicitly and not through a base class + if cfg_ana.name == "treeProducerSusyGenLevelOnly": + self.initDone = True + self.declareVariables() diff --git a/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyMultilepton.py b/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyMultilepton.py index 6870dabf517e9..b8d9221e63364 100644 --- a/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyMultilepton.py +++ b/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyMultilepton.py @@ -19,9 +19,9 @@ def __init__(self, cfg_ana, cfg_comp, looperName): NTupleVariable("htJet40a", lambda ev : ev.htJet40a, help="H_{T} computed from leptons and jets (with |eta|<4.7, pt > 40 GeV)"), NTupleVariable("mhtJet40", lambda ev : ev.mhtJet40, help="H_{T}^{miss} computed from leptons and jets (with |eta|<2.4, pt > 40 GeV)"), NTupleVariable("mhtJet40a", lambda ev : ev.mhtJet40a, help="H_{T}^{miss} computed from leptons and jets (with |eta|<4.7, pt > 40 GeV)"), - NTupleVariable("nSoftBJetLoose25", lambda ev: sum([(sv.mva>0.3 and (sv.jet == None or sv.jet.pt() < 25 or not sv.jet.btagWP("CSVM"))) for sv in ev.ivf]) + len(ev.bjetsMedium), int, help="Exclusive sum of jets with pt > 25 passing CSV medium and SV from ivf with loose sv mva"), - NTupleVariable("nSoftBJetMedium25", lambda ev: sum([(sv.mva>0.7 and (sv.jet == None or sv.jet.pt() < 25 or not sv.jet.btagWP("CSVM"))) for sv in ev.ivf]) + len(ev.bjetsMedium), int, help="Exclusive sum of jets with pt > 25 passing CSV medium and SV from ivf with medium sv mva"), - NTupleVariable("nSoftBJetTight25", lambda ev: sum([(sv.mva>0.9 and (sv.jet == None or sv.jet.pt() < 25 or not sv.jet.btagWP("CSVM"))) for sv in ev.ivf]) + len(ev.bjetsMedium), int, help="Exclusive sum of jets with pt > 25 passing CSV medium and SV from ivf with tight sv mva"), + NTupleVariable("nSoftBJetLoose25", lambda ev: sum([(sv.mva>0.3 and (sv.jet == None or sv.jet.pt() < 25)) for sv in ev.ivf]) + len(ev.bjetsMedium), int, help="Exclusive sum of jets with pt > 25 passing CSV medium and SV from ivf with loose sv mva"), + NTupleVariable("nSoftBJetMedium25", lambda ev: sum([(sv.mva>0.7 and (sv.jet == None or sv.jet.pt() < 25)) for sv in ev.ivf]) + len(ev.bjetsMedium), int, help="Exclusive sum of jets with pt > 25 passing CSV medium and SV from ivf with medium sv mva"), + NTupleVariable("nSoftBJetTight25", lambda ev: sum([(sv.mva>0.9 and (sv.jet == None or sv.jet.pt() < 25)) for sv in ev.ivf]) + len(ev.bjetsMedium), int, help="Exclusive sum of jets with pt > 25 passing CSV medium and SV from ivf with tight sv mva"), ##-------------------------------------------------- NTupleVariable("mZ1", lambda ev : ev.bestZ1[0], help="Best m(ll) SF/OS"), NTupleVariable("mZ1SFSS", lambda ev : ev.bestZ1sfss[0], help="Best m(ll) SF/SS"), @@ -114,16 +114,16 @@ def __init__(self, cfg_ana, cfg_comp, looperName): self.collections.update({ # put more here - "gentopquarks" : NTupleCollection("GenTop", genParticleType, 2, help="Generated top quarks from hard scattering"), - "genbquarks" : NTupleCollection("GenBQuark", genParticleType, 2, help="Generated bottom quarks from top quark decays"), + "gentopquarks" : NTupleCollection("GenTop", genParticleType, 6, help="Generated top quarks from hard scattering"), + "genbquarks" : NTupleCollection("GenBQuark", genParticleType, 6, help="Generated bottom quarks from top quark decays"), "genwzquarks" : NTupleCollection("GenQuark", genParticleWithSourceType, 6, help="Generated quarks from W/Z decays"), ##-------------------------------------------------- "selectedTaus" : NTupleCollection("TauGood", tauTypeSusy, 3, help="Taus after the preselection"), "otherLeptons" : NTupleCollection("LepOther", leptonTypeSusyExtra, 8, help="Leptons after the preselection"), "selectedLeptons" : NTupleCollection("LepGood", leptonTypeSusyExtra, 8, help="Leptons after the preselection"), ##------------------------------------------------ - "cleanJets" : NTupleCollection("Jet", jetTypeSusy, 8, help="Cental jets after full selection and cleaning, sorted by pt"), - "cleanJetsFwd" : NTupleCollection("JetFwd", jetTypeSusy, 4, help="Forward jets after full selection and cleaning, sorted by pt"), + "cleanJets" : NTupleCollection("Jet", jetTypeSusy, 15, help="Cental jets after full selection and cleaning, sorted by pt"), + "cleanJetsFwd" : NTupleCollection("JetFwd", jetTypeSusy, 6, help="Forward jets after full selection and cleaning, sorted by pt"), ##------------------------------------------------ "ivf" : NTupleCollection("SV", svType, 20, help="SVs from IVF"), "genBHadrons" : NTupleCollection("GenBHad", heavyFlavourHadronType, 20, mcOnly=True, help="Gen-level B hadrons"), diff --git a/CMGTools/TTHAnalysis/python/analyzers/ttHGenLevelOnlyStudy.py b/CMGTools/TTHAnalysis/python/analyzers/ttHGenLevelOnlyStudy.py new file mode 100644 index 0000000000000..185f5de9d62db --- /dev/null +++ b/CMGTools/TTHAnalysis/python/analyzers/ttHGenLevelOnlyStudy.py @@ -0,0 +1,274 @@ +import operator +import itertools +import copy + +from ROOT import TLorentzVector + +from CMGTools.RootTools.fwlite.Analyzer import Analyzer +from CMGTools.RootTools.fwlite.Event import Event +from CMGTools.RootTools.statistics.Counter import Counter, Counters +from CMGTools.RootTools.fwlite.AutoHandle import AutoHandle +from CMGTools.RootTools.physicsobjects.Lepton import Lepton +from CMGTools.RootTools.physicsobjects.Photon import Photon +from CMGTools.RootTools.physicsobjects.Electron import Electron +from CMGTools.RootTools.physicsobjects.Muon import Muon +from CMGTools.RootTools.physicsobjects.Jet import Jet +from CMGTools.RootTools.physicsobjects.PhysicsObjects import GenParticle + +from CMGTools.RootTools.utils.DeltaR import deltaR,deltaPhi +from CMGTools.RootTools.physicsobjects.genutils import * + +class LeptonFromGen: + def __init__(self, physObj): + self.physObj = physObj + def __getattr__(self, attr): + if hasattr(self.physObj, attr): + return getattr(self.physObj, attr) + raise RuntimeError, "Missing attribute '%s'" % attr +class JetFromGen: + def __init__(self, physObj): + self.physObj = physObj + self.btag = False + def btagCSV(self,tag): + return self.btag + def __getattr__(self, attr): + if hasattr(self.physObj, attr): + return getattr(self.physObj, attr) + raise RuntimeError, "Missing attribute '%s'" % attr + +class ttHGenLevelOnlyStudy( Analyzer ): + """ + Fakes a reco event starting from GEN-only files + """ + def __init__(self, cfg_ana, cfg_comp, looperName ): + super(ttHGenLevelOnlyStudy,self).__init__(cfg_ana,cfg_comp,looperName) + self.doPDFWeights = hasattr(self.cfg_ana, "PDFWeights") and len(self.cfg_ana.PDFWeights) > 0 + if self.doPDFWeights: + self.pdfWeightInit = False + #--------------------------------------------- + # DECLARATION OF HANDLES OF GEN LEVEL OBJECTS + #--------------------------------------------- + + + def declareHandles(self): + super(ttHGenLevelOnlyStudy, self).declareHandles() + self.mchandles['genParticles'] = AutoHandle( 'genParticles', 'std::vector' ) + self.mchandles['jets'] = AutoHandle( 'ak5GenJets', 'std::vector' ) + self.mchandles['met'] = AutoHandle( 'genMetTrue', 'std::vector' ) + if self.doPDFWeights: + self.mchandles['pdfstuff'] = AutoHandle( 'generator', 'GenEventInfoProduct' ) + + def beginLoop(self): + super(ttHGenLevelOnlyStudy,self).beginLoop() + + def doLeptons(self,iEvent,event): + def isPrompt(l): + for x in xrange(l.numberOfMothers()): + mom = l.mother(x) + if mom.status() > 2: return True + id = abs(mom.pdgId()) + if id > 100: return False + if id < 6: return False + if id == 21: return False + if id in [11,13,15]: return isPrompt(mom) + if id in [ 22,23,24,25,32 ]: return True + return True + + event.selectedLeptons = [] + for l in event.genParticles: + if abs(l.pdgId()) not in [11,13] or l.status() != 1: continue + if abs(l.pdgId()) == 13: + if l.pt() <= 5 or abs(l.eta()) > 2.4: continue + if abs(l.pdgId()) == 11: + if l.pt() <= 7 or abs(l.eta()) > 2.5: continue + if not isPrompt(l): + continue + event.selectedLeptons.append(LeptonFromGen(l)) + + def doJets(self,iEvent,event): + event.cleanJetsAll = [] + event.cleanJetsFwd = [] + event.cleanJets = [] + for j in self.mchandles['jets'].product(): + if j.pt() < 25: continue + for l in event.selectedLeptons: + if l.pt() > 10 and deltaR(l.eta(),l.phi(),j.eta(),j.phi()) < 0.5: + continue + jo = JetFromGen(j) + event.cleanJetsAll.append(jo) + if abs(j.eta()) < 2.4: + event.cleanJetsAll.append(jo) + else: + event.cleanJetsFwd.append(jo) + def doBTag(self,iEvent,event): + bs = [] + for gp in event.genParticles: + if gp.status() != 2: continue + id = abs(gp.pdgId()) + if id == 5 or ((id % 1000) / 100) == 5 or ((id % 10000)/1000) == 5: + bs.append(gp) + for j in event.cleanJetsAll: + for gp in bs: + if deltaR(gp.eta(),gp.phi(),j.eta(),j.phi()) < 0.4: + gp.btag = True + def doMET(self,iEvent,event): + event.met = self.mchandles['met'].product().front() + + def makeZs(self, event, maxLeps): + event.bestZ1 = [ 0., -1,-1 ] + event.bestZ1sfss = [ 0., -1,-1 ] + event.bestZ2 = [ 0., -1,-1, 0. ] + nlep = len(event.selectedLeptons) + for i,l1 in enumerate(event.selectedLeptons): + for j in range(i+1,nlep): + if j >= maxLeps: break + l2 = event.selectedLeptons[j] + if l1.pdgId() == -l2.pdgId(): + zmass = (l1.p4() + l2.p4()).M() + if event.bestZ1[0] == 0 or abs(zmass - 91.188) < abs(event.bestZ1[0] - 91.188): + event.bestZ1 = [ zmass, i, j ] + if l1.pdgId() == l2.pdgId(): + zmass = (l1.p4() + l2.p4()).M() + if event.bestZ1sfss[0] == 0 or abs(zmass - 91.188) < abs(event.bestZ1sfss[0] - 91.188): + event.bestZ1sfss = [ zmass, i, j ] + if event.bestZ1[0] != 0 and nlep > 3: + for i,l1 in enumerate(event.selectedLeptons): + if i == event.bestZ1[1]: continue + for j in range(i+1,nlep): + if j >= maxLeps: break + if j == event.bestZ1[2]: continue + l2 = event.selectedLeptons[j] + if l1.pdgId() == -l2.pdgId(): + if l1.pt() + l2.pt() > event.bestZ2[0]: + event.bestZ2 = [ l1.pt() + l2.pt(), i, j, (l1.p4() + l2.p4()).M() ] + + def makeMlls(self, event, maxLeps): + mllsfos = self.mllValues(event, lambda l1,l2 : l1.pdgId() == -l2.pdgId(), maxLeps) + mllafos = self.mllValues(event, lambda l1,l2 : l1.charge() == -l2.charge(), maxLeps) + mllafss = self.mllValues(event, lambda l1,l2 : l1.charge() == l2.charge(), maxLeps) + mllafas = self.mllValues(event, lambda l1,l2 : True, maxLeps) + event.minMllSFOS = min(mllsfos) + event.minMllAFOS = min(mllafos) + event.minMllAFSS = min(mllafss) + event.minMllAFAS = min(mllafas) + event.maxMllSFOS = max(mllsfos) + event.maxMllAFAS = max(mllafas) + event.maxMllAFOS = max(mllafos) + event.maxMllAFSS = max(mllafss) + drllafos = self.drllValues(event, lambda l1,l2 : l1.charge() == -l2.charge(), maxLeps) + drllafss = self.drllValues(event, lambda l1,l2 : l1.charge() == l2.charge(), maxLeps) + event.minDrllAFSS = min(drllafss) + event.minDrllAFOS = min(drllafos) + event.maxDrllAFOS = max(drllafos) + event.maxDrllAFSS = max(drllafss) + ptllafos = self.ptllValues(event, lambda l1,l2 : l1.charge() == -l2.charge(), maxLeps) + ptllafss = self.ptllValues(event, lambda l1,l2 : l1.charge() == l2.charge(), maxLeps) + event.minPtllAFSS = min(ptllafss) + event.minPtllAFOS = min(ptllafos) + event.maxPtllAFOS = max(ptllafos) + event.maxPtllAFSS = max(ptllafss) + leps = event.selectedLeptons; nlep = len(leps) + event.m2l = (leps[0].p4() + leps[1].p4()).M() if nlep >= 2 else 0 + event.pt2l = (leps[0].p4() + leps[1].p4()).Pt() if nlep >= 2 else 0 + event.q3l = sum([l.charge() for l in leps[:2]]) if nlep >= 3 else 0 + event.ht3l = sum([l.pt() for l in leps[:2]]) if nlep >= 3 else 0 + event.pt3l = (leps[0].p4() + leps[1].p4() + leps[2].p4()).Pt() if nlep >= 3 else 0 + event.m3l = (leps[0].p4() + leps[1].p4() + leps[2].p4()).M() if nlep >= 3 else 0 + event.q4l = sum([l.charge() for l in leps[:3]]) if nlep >= 4 else 0 + event.ht4l = sum([l.pt() for l in leps[:3]]) if nlep >= 4 else 0 + event.pt4l = (leps[0].p4() + leps[1].p4() + leps[2].p4() + leps[3].p4()).Pt() if nlep >= 4 else 0 + event.m4l = (leps[0].p4() + leps[1].p4() + leps[2].p4() + leps[3].p4()).M() if nlep >= 4 else 0 + + def mllValues(self, event, pairSelection, maxLeps): + return self.llValues(event, lambda l1,l2: (l1.p4() + l2.p4()).M(), pairSelection, maxLeps) + + def drllValues(self, event, pairSelection, maxLeps): + return self.llValues(event, lambda l1,l2: deltaR(l1.eta(), l1.phi(), l2.eta(), l2.phi()), pairSelection, maxLeps) + + def ptllValues(self, event, pairSelection, maxLeps): + return self.llValues(event, lambda l1,l2: (l1.p4() + l2.p4()).Pt(), pairSelection, maxLeps) + + def llValues(self, event, function, pairSelection, maxLeps): + pairs = [] + nlep = len(event.selectedLeptons) + for i,l1 in enumerate(event.selectedLeptons): + for j in range(i+1,nlep): + if j >= maxLeps: break + l2 = event.selectedLeptons[j] + if pairSelection(l1,l2): + pairs.append( function(l1, l2) ) + if pairs == []: pairs.append(-1) + return pairs + + + def initPDFWeights(self): + from ROOT import PdfWeightProducerTool + self.pdfWeightInit = True + self.pdfWeightTool = PdfWeightProducerTool() + for pdf in self.cfg_ana.PDFWeights: + self.pdfWeightTool.addPdfSet(pdf+".LHgrid") + self.pdfWeightTool.beginJob() + + def makePDFWeights(self, event): + if not self.pdfWeightInit: self.initPDFWeights() + self.pdfWeightTool.processEvent(self.mchandles['pdfstuff'].product()) + event.pdfWeights = {} + for pdf in self.cfg_ana.PDFWeights: + ws = self.pdfWeightTool.getWeights(pdf+".LHgrid") + event.pdfWeights[pdf] = [w for w in ws] + #print "Produced %d weights for %s: %s" % (len(ws),pdf,event.pdfWeights[pdf]) + + def doHT(self, iEvent, event): + import ROOT + + objects25 = [ j for j in event.cleanJets if j.pt() > 25 ] + event.selectedLeptons + objects30 = [ j for j in event.cleanJets if j.pt() > 30 ] + event.selectedLeptons + objects40 = [ j for j in event.cleanJets if j.pt() > 40 ] + event.selectedLeptons + objects40j = [ j for j in event.cleanJets if j.pt() > 40 ] + + event.htJet25 = sum([x.pt() for x in objects25]) + event.mhtJet25vec = ROOT.reco.Particle.LorentzVector(-1.*(sum([x.px() for x in objects25])) , -1.*(sum([x.py() for x in objects25])), 0, 0 ) + event.mhtPhiJet25 = event.mhtJet25vec.phi() + event.mhtJet25 = event.mhtJet25vec.pt() + + event.htJet30 = sum([x.pt() for x in objects30]) + event.mhtJet30vec = ROOT.reco.Particle.LorentzVector(-1.*(sum([x.px() for x in objects30])) , -1.*(sum([x.py() for x in objects30])), 0, 0 ) + event.mhtJet30 = event.mhtJet30vec.pt() + event.mhtPhiJet30 = event.mhtJet30vec.phi() + + event.htJet40 = sum([x.pt() for x in objects40]) + event.mhtJet40vec = ROOT.reco.Particle.LorentzVector(-1.*(sum([x.px() for x in objects40])) , -1.*(sum([x.py() for x in objects40])), 0, 0 ) + event.mhtJet40 = event.mhtJet40vec.pt() + event.mhtPhiJet40 = event.mhtJet40vec.phi() + + event.htJet40j = sum([x.pt() for x in objects40j]) + event.mhtJet40jvec = ROOT.reco.Particle.LorentzVector(-1.*(sum([x.px() for x in objects40j])) , -1.*(sum([x.py() for x in objects40j])), 0, 0 ) + event.mhtJet40j = event.mhtJet40jvec.pt() + event.mhtPhiJet40j = event.mhtJet40jvec.phi() + + def process(self, iEvent, event): + self.readCollections( iEvent ) + + # if not MC, nothing to do + if not self.cfg_comp.isMC: + return True + + event.genParticles = [ gp for gp in self.mchandles['genParticles'].product() ] + + event.eventWeigth = 1.0 + event.run = iEvent.eventAuxiliary().id().run() + event.lumi = iEvent.eventAuxiliary().id().luminosityBlock() + event.eventId = iEvent.eventAuxiliary().id().event() + + self.doLeptons(iEvent,event) + self.makeZs(event, 4) + self.makeMlls(event, 4) + self.doJets(iEvent,event) + self.doBTag(iEvent,event) + self.doMET(iEvent,event) + self.doHT(iEvent,event) + + # do PDF weights, if requested + if self.doPDFWeights: + self.makePDFWeights(event) + return True diff --git a/CMGTools/TTHAnalysis/python/analyzers/ttHLepTreeProducerNew.py b/CMGTools/TTHAnalysis/python/analyzers/ttHLepTreeProducerNew.py index 6c900a9af7bed..0c5642e09b280 100644 --- a/CMGTools/TTHAnalysis/python/analyzers/ttHLepTreeProducerNew.py +++ b/CMGTools/TTHAnalysis/python/analyzers/ttHLepTreeProducerNew.py @@ -26,8 +26,9 @@ def __init__(self, cfg_ana, cfg_comp, looperName): def declareHandles(self): super(ttHLepTreeProducerNew, self).declareHandles() - self.handles['TriggerResults'] = AutoHandle( ('TriggerResults','','HLT'), 'edm::TriggerResults' ) self.mchandles['GenInfo'] = AutoHandle( ('generator','',''), 'GenEventInfoProduct' ) + if hasattr(self.cfg_ana, 'triggerBits'): + self.handles['TriggerResults'] = AutoHandle( ('TriggerResults','','HLT'), 'edm::TriggerResults' ) for k,v in self.collections.iteritems(): if type(v) == tuple and isinstance(v[0], AutoHandle): self.handles[k] = v[0] @@ -90,13 +91,14 @@ def fillCoreVariables(self, tr, iEvent, event, isMC): tr.fill('evt', event.eventId) tr.fill('isData', 0 if isMC else 1) - triggerResults = self.handles['TriggerResults'].product() - for T,TC in self.triggerBitCheckers: - tr.fill("HLT_"+T, TC.check(iEvent.object(), triggerResults)) + if self.triggerBitCheckers: + triggerResults = self.handles['TriggerResults'].product() + for T,TC in self.triggerBitCheckers: + tr.fill("HLT_"+T, TC.check(iEvent.object(), triggerResults)) if isMC: ## PU weights - tr.fill("nTrueInt", event.nPU) + tr.fill("nTrueInt", getattr(event, 'nPU', -1)) tr.fill("puWeight", event.eventWeight) tr.fill("genWeight", self.mchandles['GenInfo'].product().weight()) ## PDF weights diff --git a/CMGTools/TTHAnalysis/python/analyzers/ttHSVAnalyzer.py b/CMGTools/TTHAnalysis/python/analyzers/ttHSVAnalyzer.py index 8cf227f8677fe..3f4b244160aa4 100644 --- a/CMGTools/TTHAnalysis/python/analyzers/ttHSVAnalyzer.py +++ b/CMGTools/TTHAnalysis/python/analyzers/ttHSVAnalyzer.py @@ -54,6 +54,7 @@ def process(self, iEvent, event): # attach distances to PV pv = event.goodVertices[0] for sv in allivf: + #sv.dz = SignedImpactParameterComputer.vertexDz(sv, pv) sv.dxy = SignedImpactParameterComputer.vertexDxy(sv, pv) sv.d3d = SignedImpactParameterComputer.vertexD3d(sv, pv) sv.cosTheta = SignedImpactParameterComputer.vertexDdotP(sv, pv) @@ -140,12 +141,24 @@ def ref2id(ref): for i in xrange(s.numberOfDaughters()): daumap[ref2id(s.daughterPtr(i))] = s for j in event.jetsIdOnly: + #print "jet with pt %5.2f, eta %+4.2f, phi %+4.2f: " % (j.pt(), j.eta(), j.phi()) jdaus = [ref2id(j.daughterPtr(i)) for i in xrange(j.numberOfDaughters())] j.svs = [] for jdau in jdaus: if jdau in daumap: + #print " --> matched by ref with SV with pt %5.2f, eta %+4.2f, phi %+4.2f: " % (daumap[jdau].pt(), daumap[jdau].eta(), daumap[jdau].phi()) j.svs.append(daumap[jdau]) daumap[jdau].jet = j + for s in event.ivf: + if s.jet != None: continue + #print "Unassociated SV with %d tracks, mass %5.2f, pt %5.2f, eta %+4.2f, phi %+4.2f: " % (s.numberOfDaughters(), s.mass(), s.pt(), s.eta(), s.phi()) + bestDr = 0.4 + for j in event.jetsIdOnly: + dr = deltaR(s.eta(),s.phi(),j.eta(),j.phi()) + if dr < bestDr: + bestDr = dr + s.jet = j + #print " close to jet with pt %5.2f, eta %+4.2f, phi %+4.2f: dr = %.3f" % (j.pt(), j.eta(), j.phi(), dr) diff --git a/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_ee.txt b/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_ee.txt index 6bce660375cb5..0ddc641c77486 100644 --- a/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_ee.txt +++ b/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_ee.txt @@ -1,7 +1,7 @@ == 2 good leptons: nLepGood10 == 2 cleanup: minMllAFAS > 12 pt2010: LepGood1_pt>20 && LepGood2_pt>20 -lep MVA: min(LepGood1_mvaNew,LepGood2_mvaNew) > 0.93 +lep MVA: min(LepGood1_mvaSusy,LepGood2_mvaSusy) > 0.93 el el: abs(LepGood1_pdgId) == 11 && abs(LepGood2_pdgId) == 11 same-sign: LepGood1_charge*LepGood2_charge > 0 tight-charge: LepGood1_tightCharge > 1 && LepGood2_tightCharge > 1 diff --git a/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_em.txt b/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_em.txt index f020442f26420..b16c0b8eb30d5 100644 --- a/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_em.txt +++ b/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_em.txt @@ -1,7 +1,7 @@ == 2 good leptons: nLepGood10 == 2 cleanup: minMllAFAS > 12 pt2010: LepGood1_pt>20 && LepGood2_pt>20 -lep MVA: min(LepGood1_mvaNew,LepGood2_mvaNew) > 0.93 +lep MVA: min(LepGood1_mvaSusy,LepGood2_mvaSusy) > 0.93 el mu: abs(LepGood1_pdgId) != abs(LepGood2_pdgId) same-sign: LepGood1_charge*LepGood2_charge > 0 tight-charge: LepGood1_tightCharge > (abs(LepGood1_pdgId) == 11) && LepGood2_tightCharge > (abs(LepGood2_pdgId) == 11) diff --git a/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_mumu.txt b/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_mumu.txt index aad52e83e836f..bf5076cf70425 100644 --- a/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_mumu.txt +++ b/CMGTools/TTHAnalysis/python/plotter/bins/susy_2lss_mumu.txt @@ -1,7 +1,7 @@ == 2 good leptons: nLepGood10 == 2 cleanup: minMllAFAS > 12 pt2010: LepGood1_pt>20 && LepGood2_pt>20 -lep MVA: min(LepGood1_mvaNew,LepGood2_mvaNew) > 0.93 +lep MVA: min(LepGood1_mvaSusy,LepGood2_mvaSusy) > 0.93 mu mu: abs(LepGood1_pdgId) == 13 && abs(LepGood2_pdgId) == 13 same-sign: LepGood1_charge*LepGood2_charge > 0 tight-charge: LepGood1_tightCharge && LepGood2_tightCharge diff --git a/CMGTools/TTHAnalysis/python/plotter/functions.cc b/CMGTools/TTHAnalysis/python/plotter/functions.cc index 3077a712863e7..4a940cffda0e3 100644 --- a/CMGTools/TTHAnalysis/python/plotter/functions.cc +++ b/CMGTools/TTHAnalysis/python/plotter/functions.cc @@ -76,6 +76,14 @@ float mt_llv(float ptl1, float phil1, float ptl2, float phil2, float ptv, float return std::sqrt(std::max(0.f, ht*ht - px*px - py*py)); } +float mt_lllv(float ptl1, float phil1, float ptl2, float phil2, float ptl3, float phil3, float ptv, float phiv) { + float px = ptl1*std::cos(phil1) + ptl2*std::cos(phil2) + ptl3*std::cos(phil3) + ptv*std::cos(phiv); + float py = ptl1*std::sin(phil1) + ptl2*std::sin(phil2) + ptl3*std::sin(phil3) + ptv*std::sin(phiv); + float ht = ptl1+ptl2+ptl3+ptv; + return std::sqrt(std::max(0.f, ht*ht - px*px - py*py)); +} + + float mtw_wz3l(float pt1, float eta1, float phi1, float m1, float pt2, float eta2, float phi2, float m2, float pt3, float eta3, float phi3, float m3, float mZ1, float met, float metphi) { if (abs(mZ1 - mass_2(pt1,eta1,phi1,m1,pt2,eta2,phi2,m2)) < 0.01) return mt_2(pt3,phi3,met,metphi); diff --git a/CMGTools/TTHAnalysis/python/plotter/mcAnalysis.py b/CMGTools/TTHAnalysis/python/plotter/mcAnalysis.py index 11b06e549570a..ac1f9fd9f50c3 100755 --- a/CMGTools/TTHAnalysis/python/plotter/mcAnalysis.py +++ b/CMGTools/TTHAnalysis/python/plotter/mcAnalysis.py @@ -273,7 +273,7 @@ def prettyPrint(self,reports,makeSummary=True): elif len(h) <= fmtlen: print h.center(fmtlen), else: - print h[:fmtlen] + print h[:fmtlen], print "" print "-"*((fmtlen+1)*len(table)+clen) for i,(cut,dummy) in enumerate(table[0][1]): diff --git a/CMGTools/TTHAnalysis/python/plotter/mcPlots.py b/CMGTools/TTHAnalysis/python/plotter/mcPlots.py index 8d24cfea436a0..8bb384591649b 100644 --- a/CMGTools/TTHAnalysis/python/plotter/mcPlots.py +++ b/CMGTools/TTHAnalysis/python/plotter/mcPlots.py @@ -477,6 +477,7 @@ def run(self,mca,cuts,plots,makeStack=True,makeCanvas=True): # # blinding policy blind = pspec.getOption('Blinded','None') if 'data' in pmap else 'None' + if options.unblind == True: blind = 'None' xblind = [9e99,-9e99] if re.match(r'(bin|x)\s*([<>]?)\s*(\+|-)?\d+(\.\d+)?|(\+|-)?\d+(\.\d+)?\s*<\s*(bin|x)\s*<\s*(\+|-)?\d+(\.\d+)?', blind): xfunc = (lambda h,b: b) if 'bin' in blind else (lambda h,b : h.GetXaxis().GetBinCenter(b)); @@ -711,6 +712,7 @@ def addPlotMakerOptions(parser): parser.add_option("--plotmode", dest="plotmode", type="string", default="stack", help="Show as stacked plot (stack), a non-stacked comparison (nostack) and a non-stacked comparison of normalized shapes (norm)") parser.add_option("--rebin", dest="globalRebin", type="int", default="0", help="Rebin all plots by this factor") parser.add_option("--poisson", dest="poisson", action="store_true", default=False, help="Draw Poisson error bars") + parser.add_option("--unblind", dest="unblind", action="store_true", default=False, help="Unblind plots irrespectively of plot file") parser.add_option("--select-plot", "--sP", dest="plotselect", action="append", default=[], help="Select only these plots out of the full file") parser.add_option("--exclude-plot", "--xP", dest="plotexclude", action="append", default=[], help="Exclude these plots from the full file") parser.add_option("--legendWidth", dest="legendWidth", type="float", default=0.25, help="Width of the legend") diff --git a/CMGTools/TTHAnalysis/python/plotter/mca-CSA14.txt b/CMGTools/TTHAnalysis/python/plotter/mca-CSA14.txt index 5749dc9f374f7..2197b7b8ad279 100644 --- a/CMGTools/TTHAnalysis/python/plotter/mca-CSA14.txt +++ b/CMGTools/TTHAnalysis/python/plotter/mca-CSA14.txt @@ -32,16 +32,16 @@ WJets : WJetsToLNu_HT400to600_PU_S14_POSTLS170 : 55.61 * 1.23; Label="W+je WJets : WJetsToLNu_HT600toInf_PU_S14_POSTLS170 : 18.81 * 1.23; Label="W+jets", FillColor=ROOT.kCyan+2, NormSystematic=0.5 # -- SUSY -T1tttt_HL+: T1tttt2J_6_PU_S14_POSTLS170 : 0.0141903 ; FillColor=ROOT.kAzure+8, Label="T1t^{4} 1.5/0.1" -T1tttt_HM+: T1tttt2J_7_PU_S14_POSTLS170 : 0.0856418 ; FillColor=ROOT.kAzure+0, Label="T1t^{4} 1.2/0.8" +T1tttt_HL+: SMS_T1tttt_2J_mGl1500_mLSP100_PU_S14_POSTLS170 : 0.0141903 ; FillColor=ROOT.kAzure+8, Label="T1t^{4} 1.5/0.1" +T1tttt_HM+: SMS_T1tttt_2J_mGl1200_mLSP800_PU_S14_POSTLS170 : 0.0856418 ; FillColor=ROOT.kAzure+0, Label="T1t^{4} 1.2/0.8" T5Full_HM+: T5Full_1200_1000_800 : 0.0856418 ; FillColor=ROOT.kViolet+6, Label="T5W^{2} 1.2/0.8" T5Full_HL+: T5Full_1500_800_100 : 0.0141903 ; FillColor=ROOT.kViolet+0, Label="T5W^{2} 1.5/1.1" -GEN_T5tttt_HLDx: GEN_T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_py8 : 0.0460525 ; FillColor=ROOT.kAzure+8, Label="G5t^{4} 1.3/.3/\#chi", TreeName="treeProducerSusyGenLevelOnly" -GEN_T5tttt_HLD4: GEN_T1tttt_2J_mGo1300_mStop300_mChi280_py8 : 0.0460525 ; FillColor=ROOT.kRed-9, Label="G5t^{4} 1.3/.3/4", TreeName="treeProducerSusyGenLevelOnly" -GEN_T5tttt_LLDx: GEN_T1tttt_2J_mGo800_mStop300_mCh285_mChi280_py8 : 1.4891 ; FillColor=ROOT.kAzure+8, Label="G5t^{4} 0.8/.3/\#chi", TreeName="treeProducerSusyGenLevelOnly" -GEN_T5tttt_LLD4: GEN_T1tttt_2J_mGo800_mStop300_mChi280_py8 : 1.4891 ; FillColor=ROOT.kRed-9, Label="G5t^{4} 0.8/.3/4", TreeName="treeProducerSusyGenLevelOnly" +#GEN_T5tttt_HLDx: GEN_T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_py8 : 0.0460525 ; FillColor=ROOT.kAzure+8, Label="G5t^{4} 1.3/.3/\#chi", TreeName="treeProducerSusyGenLevelOnly" +#GEN_T5tttt_HLD4: GEN_T1tttt_2J_mGo1300_mStop300_mChi280_py8 : 0.0460525 ; FillColor=ROOT.kRed-9, Label="G5t^{4} 1.3/.3/4", TreeName="treeProducerSusyGenLevelOnly" +#GEN_T5tttt_LLDx: GEN_T1tttt_2J_mGo800_mStop300_mCh285_mChi280_py8 : 1.4891 ; FillColor=ROOT.kAzure+8, Label="G5t^{4} 0.8/.3/\#chi", TreeName="treeProducerSusyGenLevelOnly" +#GEN_T5tttt_LLD4: GEN_T1tttt_2J_mGo800_mStop300_mChi280_py8 : 1.4891 ; FillColor=ROOT.kRed-9, Label="G5t^{4} 0.8/.3/4", TreeName="treeProducerSusyGenLevelOnly" T5tttt_HLDx+: T1tttt_2J_mGo1300_mStop300_mCh285_mChi280_pythia8_S14: 0.0460525 ; FillColor=ROOT.kRed-9, Label="T5t^{4} 1.3/.3/\#chi" T5tttt_HLD4+: T1tttt_2J_mGo1300_mStop300_mChi280_pythia8_S14 : 0.0460525 ; FillColor=ROOT.kRed+0, Label="T5t^{4} 1.3/.3/4" T5tttt_LLDx+: T1tttt_2J_mGo800_mStop300_mCh285_mChi280_pythia8_S14 : 1.4891 ; FillColor=ROOT.kMagenta-9, Label="T5t^{4} 0.8/.3/\#chi" diff --git a/CMGTools/TTHAnalysis/python/plotter/simpleOptimizer.py b/CMGTools/TTHAnalysis/python/plotter/simpleOptimizer.py index c864a53052694..ded6072b96246 100644 --- a/CMGTools/TTHAnalysis/python/plotter/simpleOptimizer.py +++ b/CMGTools/TTHAnalysis/python/plotter/simpleOptimizer.py @@ -20,6 +20,7 @@ def getPoints(self): class ParamSampler: def __init__(self,params): self._params = params + self._cache = {} def _addAllPoints(self, param, points): ## trivial case if len(points) == 0: @@ -38,21 +39,180 @@ def getAllPoints(self): return [ dict(x) for x in allpoints ] def optimize(self,mca,cut,fom,algo): if algo == "scan": return self.scan(mca,cut,fom,algo) + if algo == "walk": return self.walk(mca,cut,fom,algo) + if "cloud:" in algo: return self.cloud(mca,cut,fom,algo) + if "scipy:" in algo: return self.scipyOpt(mca,cut,fom,algo) + def onepoint(self,mca,cut,fom,point): + if self._cache != None: + key = repr(point) + if key in self._cache: + return self._cache[key] + mycut = CutsFile([('all',cut.format(**point))]) ## need the '**' to unpack a map + yields = mca.getYields(mycut,nodata=True,noEntryLine=True) + ret = fom(yields) + if self._cache != None: + self._cache[repr(point)] = ret + return ret def scan(self,mca,cut,fom,algo): best, bval = None, None allpoints = self.getAllPoints() print "I have %d points to test... please sit back and relax." % len(allpoints) for i,point in enumerate(allpoints): - print "trying point ",i,": ",point - mycut = CutsFile([('all',cut.format(**point))]) ## need the '**' to unpack a map - yields = mca.getYields(mycut,nodata=True,noEntryLine=True) - fomval = fom(yields) + print "trying point ",i,": ",point, + fomval = self.onepoint(mca,cut,fom,point) print " ---> ",fomval - if best == None or bval < fomval: + if best == None or bval[0] < fomval[0]: best = point bval = fomval return (best,bval) - + def scan1(self,mca,cut,fom,algo,param,cursor,value): + print "Scanning %s values [ %s ] starting from %s" % (param.name, param.getPoints(), cursor ) + best, bval = cursor[param.name], value + for pv in param.getPoints(): + if pv == best: continue + cursor[param.name] = pv + fomval = self.onepoint(mca,cut,fom,cursor) + if fomval[0] > bval[0]: + best = pv; bval = fomval + cursor[param.name] = best; + return (cursor, best) + def expand(self,mca,cut,fom,cursor,value,size=1): + neighbours = [ [] ] # list of points, each of which is a lists of pairs (name,value) + for param in self._params: + values = param.getPoints() + index = values.index(cursor[param.name]) + newneigh = [] + for i in xrange(max(0,index-size),min(index+size,len(values)-1)+1): + for prefix in neighbours: + newneigh.append( prefix + [ (param.name, values[i]) ] ) + neighbours = newneigh + neighbours = [ dict(x) for x in neighbours ] + print "Scanning all the %d neighbours within %d: " % (len(neighbours), size) + best, bval = cursor, value + for i,point in enumerate(neighbours): + fomval = self.onepoint(mca,cut,fom,point) + if fomval[0] > bval[0]: + print " %6d ---> %s (improvement)" % (i,fomval) + best = point; bval = fomval + if bval[0] > value[0]: + print " found new best point %s fom value %s" % (best, bval) + return (best, bval, bval[0] > value[0]) + def step1(self,mca,cut,fom,algo,param,cursor,value): + print "Stepping %s values [ %s ] starting from %s at %s" % (param.name, param.getPoints(), cursor, value ) + best, bval = cursor[param.name], value + values = param.getPoints() + index = values.index(cursor[param.name]) + ileft = index - 1; moved = False + while ileft >= 0: + cursor[param.name] = values[ileft] + fomval = self.onepoint(mca,cut,fom,cursor) + if fomval[0] > bval[0]: + best = values[ileft]; bval = fomval; moved = True + ileft -= 1 + else: + break + if moved: + cursor[param.name] = best; + print " moved %s to %s, fom value %s" % (param.name, best, bval) + return (cursor, bval, True) + iright = index + 1 + while iright < len(values): + cursor[param.name] = values[iright] + fomval = self.onepoint(mca,cut,fom,cursor) + if fomval[0] > bval[0]: + best = values[iright]; bval = fomval; moved = True + iright += 1 + else: + break + if moved: + cursor[param.name] = best; + print " moved %s to %s, fom value %s" % (param.name, best, bval) + return (cursor, bval, True) + cursor[param.name] = best; + print " did not move %s from %s at %s" % (param.name, best, bval) + return (cursor, bval, False) + def walk(self,mca,cut,fom,algo): + cursor = dict( [ (p.name,p.value) for p in self._params ] ) + fomval = self.onepoint(mca,cut,fom,cursor) + print "Starting point: %s (fomval = %s)" % (cursor, fomval) + active = dict([ (p.name,True) for p in self._params ]) + constant = dict([ (p.name,len(p.getPoints()) == 1) for p in self._params ]) + for iter in xrange(1000): + walked = False + lastactive = {} + for p in self._params: + if constant[p.name]: continue + if active[p.name]: + lastactive[p.name] = True + (cursor, fomval, act) = self.step1(mca,cut,fom,algo,p,cursor,fomval) + active[p.name] = act + if act: walked = True + if walked: continue + for p in self._params: + if constant[p.name]: continue + if p.name in lastactive: continue # these are already at the maximum + (cursor, fomval, act) = self.step1(mca,cut,fom,algo,p,cursor,fomval) + active[p.name] = act + if act: walked = True + if walked: + for pn in active.keys(): active[pn] = True + else: + break + return cursor, fomval + def cloud(self,mca,cut,fom,algo): + cursor = dict( [ (p.name,p.value) for p in self._params ] ) + fomval = self.onepoint(mca,cut,fom,cursor) + print "Starting point: %s (fomval = %s)" % (cursor, fomval) + for iter in xrange(1000): + (cursor, fomval, moved) = self.expand(mca,cut,fom,cursor,fomval,size=int(algo.replace("cloud:",""))) + if not moved: + break + return (cursor,fomval) + def scipyOpt(self,mca,cut,fom,algo): + cursor = dict( [ (p.name,p.value) for p in self._params ] ) + fomval = self.onepoint(mca,cut,fom,cursor) + print "Starting point: %s (fomval = %s)" % (cursor, fomval) + import scipy.optimize + minimizer = getattr(scipy.optimize,algo.replace("scipy:","")) + fun = ParamSampelFunctor(self,mca,cut,fom) + ret = minimizer(fun, fun.x0) + cursor = fun.int2ext(ret) + fomval = self.onepoint(mca,cut,fom,cursor) + print "Ending point: %s (fomval = %s)" % (cursor, fomval) + return (cursor,fomval) +class ParamSampelFunctor: + def __init__(self,paramSampler,mca,cut,fom): + self._sampler = paramSampler + self._mca = mca + self._cut = cut + self._fom = fom + self._i = 0 + self.xnames = [] + self.xmin = [] + self.xmax = [] + self.x0 = [] + self.consts = [] + for p in self._sampler._params: + points = p.getPoints() + if len(points) > 1: + self.xnames.append(p.name) + self.xmax.append(max(points)) + self.xmin.append(min(points)) + self.x0.append(0.5) + else: + self.consts.append((p.name,p.value)) + def int2ext(self,x): + point = dict(self.consts) + for (n,xmin,xmax,xraw) in zip(self.xnames,self.xmin,self.xmax,x): + point[n] = max(0,0.5*((xraw-0.5)*10+1))*(xmax-xmin) + xmin + return point + def __call__(self,x): + point = self.int2ext(x) + fomval = self._sampler.onepoint(self._mca,self._cut,self._fom, point) + self._i += 1 + print " %6d: %s --> %s " % (self._i, point,fomval) + return -fomval[0] + class ParamFile: def __init__(self,params,options): self._params = [] @@ -91,8 +251,22 @@ def __call__(self,report): totb += report[b][-1][1][0] totbsyst += (syst*report[b][-1][1][0])**2 if self._includeS: totb += tots - return tots/sqrt(totb + totbsyst**2) - + return (tots/sqrt(totb + totbsyst**2), "S=%.2f, B=%.2f" % (tots, totb)) +class PunziFOM: + def __init__(self,mca,syst=0.0,a=3): + self._signals = mca.listSignals() + self._backgrounds = [ (b,syst) for b in mca.listBackgrounds() ] + self._a = a + def __call__(self,report): + tots, totb, totbsyst = 0.,0.,0. + for s in self._signals: + if s not in report: continue + tots += report[s][-1][1][0] + for b,syst in self._backgrounds: + if b not in report: continue + totb += report[b][-1][1][0] + totbsyst += (syst*report[b][-1][1][0])**2 + return (tots/(0.5*self._a + sqrt(totb + totbsyst**2)), "S=%.2f, B=%.2f" % (tots, totb)) if __name__ == "__main__": from optparse import OptionParser @@ -101,7 +275,8 @@ def __call__(self,report): (options, args) = parser.parse_args() mca = MCAnalysis(args[0],options) cuts = CutsFile(args[1],options) - fom = SimpleFOM(mca) + #fom = SimpleFOM(mca) + fom = PunziFOM(mca,a=2) algo = args[2] params = ParamFile(args[3],options) sampler = ParamSampler(params.allParams()) diff --git a/CMGTools/TTHAnalysis/python/samples/ComponentCreator.py b/CMGTools/TTHAnalysis/python/samples/ComponentCreator.py index 3b4277c667d08..4de0626d785bb 100644 --- a/CMGTools/TTHAnalysis/python/samples/ComponentCreator.py +++ b/CMGTools/TTHAnalysis/python/samples/ComponentCreator.py @@ -5,12 +5,12 @@ import re class ComponentCreator(object): - def makeMCComponent(self,name,dataset,user,pattern): + def makeMCComponent(self,name,dataset,user,pattern,useAAA=False): component = cfg.MCComponent( dataset=dataset, name = name, - files = self.getFiles(dataset,user,pattern), + files = self.getFiles(dataset,user,pattern,useAAA=useAAA), xSection = 1, nGenEvents = 1, triggers = [], @@ -41,12 +41,12 @@ def makePrivateMCComponent(self,name,dataset,files): return component ### MM - def makeMyPrivateMCComponent(self,name,dataset,user,pattern,dbsInstance): + def makeMyPrivateMCComponent(self,name,dataset,user,pattern,dbsInstance, useAAA=False): component = cfg.MCComponent( dataset=dataset, name = name, - files = self.getMyFiles(dataset, user, pattern, dbsInstance), + files = self.getMyFiles(dataset, user, pattern, dbsInstance, useAAA=useAAA), xSection = 1, nGenEvents = 1, triggers = [], @@ -56,7 +56,7 @@ def makeMyPrivateMCComponent(self,name,dataset,user,pattern,dbsInstance): return component ### MM - def makeMCComponentFromEOS(self,name,dataset,path,pattern=".*root"): + def getFilesFromEOS(self,name,dataset,path,pattern=".*root"): from CMGTools.Production.dataset import getDatasetFromCache, writeDatasetToCache if "%" in path: path = path % dataset; try: @@ -66,10 +66,12 @@ def makeMCComponentFromEOS(self,name,dataset,path,pattern=".*root"): if len(files) == 0: raise RuntimeError, "ERROR making component %s: no files found under %s matching '%s'" % (name,path,pattern) writeDatasetToCache('EOS%{path}%{pattern}.pck'.format(path = path.replace('/','_'), pattern = pattern), files) + return files + def makeMCComponentFromEOS(self,name,dataset,path,pattern=".*root"): component = cfg.MCComponent( dataset=dataset, name = name, - files = files, + files = self.getFilesFromEOS(name,dataset,path,pattern), xSection = 1, nGenEvents = 1, triggers = [], @@ -97,18 +99,22 @@ def makeDataComponent(self,name,datasets,user,pattern): return component - def getFiles(self, dataset, user, pattern): + def getFiles(self, dataset, user, pattern, useAAA=False): # print 'getting files for', dataset,user,pattern ds = datasetToSource( user, dataset, pattern, True ) files = ds.fileNames - return ['root://eoscms.cern.ch//eos/cms%s' % f for f in files] + mapping = 'root://eoscms.cern.ch//eos/cms%s' + if useAAA: mapping = 'root://cms-xrd-global.cern.ch/%s' + return [ mapping % f for f in files] - def getMyFiles(self, dataset, user, pattern, dbsInstance): + def getMyFiles(self, dataset, user, pattern, dbsInstance, useAAA=False): # print 'getting files for', dataset,user,pattern ds = myDatasetToSource( user, dataset, pattern, dbsInstance, True ) files = ds.fileNames - return ['root://eoscms.cern.ch//eos/cms%s' % f for f in files] + mapping = 'root://eoscms.cern.ch//eos/cms%s' + if useAAA: mapping = 'root://cms-xrd-global.cern.ch/%s' + return [ mapping % f for f in files] def getSkimEfficiency(self,dataset,user): info=DatasetInformation(dataset,user,'',False,False,'','','') diff --git a/CMGTools/TTHAnalysis/python/samples/samples_13TeV_CSA14.py b/CMGTools/TTHAnalysis/python/samples/samples_13TeV_CSA14.py index b9d323111b8eb..7c7dc290c04a4 100644 --- a/CMGTools/TTHAnalysis/python/samples/samples_13TeV_CSA14.py +++ b/CMGTools/TTHAnalysis/python/samples/samples_13TeV_CSA14.py @@ -323,7 +323,10 @@ SingleMu = cfg.DataComponent( name = 'SingleMu', - files = [ 'root://eoscms//eos/cms/store/cmst3/user/gpetrucc/miniAOD/v1/SingleMu-Run2012D-15Apr2014-v1_PAT_big.root' ], + files = kreator.getFilesFromEOS('SingleMu', + '/SingleMu/Run2012D-15Apr2014-v1/AOD/02e0a1be-c9c7-11e3-bfe2-0024e83ef644/MINIAOD/CMSSW_7_0_9_patch2_GR_70_V2_AN1', + '/eos/cms/store/cmst3/user/cmgtools/CMG/%s'), + #[ 'root://eoscms//eos/cms/store/cmst3/user/gpetrucc/miniAOD/v1/SingleMu-Run2012D-15Apr2014-v1_PAT_big.root' ], intLumi = 1, triggers = [], json = json @@ -357,6 +360,6 @@ comp.efficiency = eff2012 for comp in dataSamplesAll: - comp.splitFactor = 1 + comp.splitFactor = 1000 comp.isMC = False comp.isData = True diff --git a/CMGTools/TTHAnalysis/python/treeReAnalyzer.py b/CMGTools/TTHAnalysis/python/treeReAnalyzer.py index 712a2eb5fd1d2..9762a1d940ee5 100644 --- a/CMGTools/TTHAnalysis/python/treeReAnalyzer.py +++ b/CMGTools/TTHAnalysis/python/treeReAnalyzer.py @@ -65,9 +65,19 @@ def __init__(self,event,prefix,index=None): def __getattr__(self,name): if name in self.__dict__: return self.__dict__[name] if name == "pdgLabel": return self.pdgLabel_() - val = getattr(self._event,self._prefix+name) - if self._index != None: - val = val[self._index] + try: + val = getattr(self._event,self._prefix+name) + if self._index != None: + val = val[self._index] + except AttributeError, e: + import math + if hasattr(math, name): val = getattr(math,name) + elif hasattr(__builtins__,name): val = getattr(__builtins__,name) + else: + try: + val = getattr(ROOT,name) + except AttributeError, e2: + raise e self.__dict__[name] = val ## cache return val def __getitem__(self,attr):