From 101d3094175a92963195ef18e7bb25248d9590ef Mon Sep 17 00:00:00 2001 From: Salvatore Rappoccio Date: Tue, 5 Nov 2013 13:53:34 -0600 Subject: [PATCH] As per request, moving all of the IHeartNY scripts to a new directory --- EDSHyFT/test/iheartny_topxs_fwlite.py | 1197 ------------------------- 1 file changed, 1197 deletions(-) delete mode 100644 EDSHyFT/test/iheartny_topxs_fwlite.py diff --git a/EDSHyFT/test/iheartny_topxs_fwlite.py b/EDSHyFT/test/iheartny_topxs_fwlite.py deleted file mode 100644 index c3d1121b..00000000 --- a/EDSHyFT/test/iheartny_topxs_fwlite.py +++ /dev/null @@ -1,1197 +0,0 @@ -#! /usr/bin/env python -import os -import glob -import time - - -from optparse import OptionParser - - -parser = OptionParser() - -parser.add_option('--files', metavar='F', type='string', action='store', - dest='files', - help='Input files') - - -parser.add_option('--outname', metavar='F', type='string', action='store', - default='TTSemilepAnalyzer_antibtag_w_mucut', - dest='outname', - help='output name') - -parser.add_option('--set', metavar='F', type='string', action='store', - default='madgraph', - dest='set', - help='mcatnlo or powheg') - -parser.add_option('--num', metavar='F', type='string', action='store', - default = 'all', - dest = 'num', - help = 'file number events in millions') - -parser.add_option('--etabin', metavar='F', type='string', action='store', - default = 'all', - dest = 'etabin', - help = 'range1 range2 or all') - - -parser.add_option('--pileup', metavar='F', type='string', action='store', - default='none', - dest='pileup', - help='ttbar or wjets') - -parser.add_option('--ptWeight', metavar='F', action='store_true', - default=False, - dest='ptWeight', - help='do pt reweighting (for MC)') - -parser.add_option('--type2', metavar='F', action='store_true', - default=False, - dest='type2', - help='look at type 2 top events') - -parser.add_option('--muOnly', metavar='F', action='store_true', - default=False, - dest='muOnly', - help='use only muons') - -parser.add_option('--useClosestForTopMass', metavar='F', action='store_true', - default=False, - dest='useClosestForTopMass', - help='use closest jet for top mass, instead of b-jet in had. hemisph.') - - -parser.add_option('--useLoose', metavar='F', action='store_true', - default=False, - dest='useLoose', - help='use loose leptons (exclusive from tight)') - - -parser.add_option('--printEvents', metavar='F', action='store_true', - default=False, - dest='printEvents', - help='Print events that pass selection (run:lumi:event)') - - -parser.add_option('--htCut', metavar='F', type='float', action='store', - default=None, - dest='htCut', - help='HT cut') - -parser.add_option('--makeResponse', metavar='M', action='store_true', - default=False, - dest='makeResponse', - help='Make response for top pt unfolding') - - -parser.add_option('--ptCut', metavar='F', type='float', action='store', - default=200.0, - dest='ptCut', - help='Leading jet PT cut') - -parser.add_option('--bDiscCut', metavar='F', type='float', action='store', - default=0.679, - dest='bDiscCut', - help='B discriminator cut') - -(options, args) = parser.parse_args() - -argv = [] - -eventsbegin = [1,10000001,20000001,30000001,40000001,50000001,60000001,70000001] -eventsend = [10000000,20000000,30000000,40000000,50000000,60000000,70000000,80000000] - -if options.num != 'all': - ifile=int(options.num) - -import ROOT -ROOT.gROOT.Macro("rootlogon.C") - - -if options.makeResponse : - ROOT.gSystem.Load("RooUnfold-1.1.1/libRooUnfold") - - -# Define classes that use ROOT - -def findClosestJet( ijet, hadJets ) : - minDR = 9999. - ret = -1 - for jjet in range(0,len(hadJets) ): - if ijet == jjet : - continue - dR = hadJets[ijet].DeltaR(hadJets[jjet]) - if dR < minDR : - minDR = hadJets[ijet].DeltaR(hadJets[jjet]) - ret = jjet - return ret - - -def findClosestInList( p41, p4list ) : - minDR = 9999. - ret = None - for j in range(0,len(p4list) ): - dR = p4list[ijet].DeltaR(p41) - if dR < minDR : - minDR = dR - ret = p4list[ijet] - return ret - - -class GenTopQuark : - pdgId = 6 # 6 = top, -6 = antitop - p4 = ROOT.TLorentzVector() - decay = 0 # 0 = hadronic, 1 = leptonic - def __init__( self, pdgId, p4, decay ) : - self.pdgId = pdgId - self.p4 = p4 - self.decay = decay - def match( self, jets ) : - return findClosestInList( self.p4, jets ) - - - - - - - - -import sys -from DataFormats.FWLite import Events, Handle - - -start_time = time.time() - - -print 'Getting these files : ' + options.files - -files = glob.glob( options.files ) -print files -toptype = '_type1' -ptwstring='' -if options.type2 == True : - toptype='_type2' -if options.ptWeight == True : - ptwstring='_ptWeighted' -etastring = '' -if options.etabin == 'range1': - etastring = '_eta_range1' -if options.etabin == 'range2': - etastring = '_eta_range2' - -if options.num != 'all': - f = ROOT.TFile( "partialfiles/"+options.outname +options.num+ptwstring+toptype+etastring+".root", "recreate" ) - name = options.outname +options.num+ptwstring+toptype+etastring -else: - f = ROOT.TFile(options.outname + ptwstring+toptype+etastring+".root", "recreate") - name = options.outname+ptwstring+toptype+etastring - - -print "Creating histograms" - -PileFile = ROOT.TFile("Pileup_plots.root") -PilePlot = PileFile.Get("pweight" + options.pileup) - -f.cd() -nJets = ROOT.TH1F("nJets", "Number of Jets, p_{T} > 30 GeV;N_{Jets};Number", 20, -0.5, 19.5 ) -nMuons = ROOT.TH1F("nMuons", "Number of Muons, p_{T} > 45 GeV;N_{Muons};Number", 5, -0.5, 4.5 ) -nElectrons = ROOT.TH1F("nElectrons", "Number of Electrons, p_{T} > 60 GeV;N_{Jets};Number", 5, -0.5, 4.5 ) - -ptMu = ROOT.TH1F("ptMu", "p_{T} of Muon", 200, 0., 200.) -ptMET0 = ROOT.TH1F("ptMET0", "MET", 200, 0., 200.) -ptMET1 = ROOT.TH1F("ptMET1", "MET", 200, 0., 200.) -ptMET2 = ROOT.TH1F("ptMET2", "MET", 200, 0., 200.) -ptMET3 = ROOT.TH1F("ptMET3", "MET", 200, 0., 200.) -htLep3mu = ROOT.TH1F("htlep3mu", "H_{T}^{Lep}", 300, 0., 600.) -htLep3t1kin = ROOT.TH1F("htlep3t1kin", "H_{T}^{Lep}", 300, 0., 600.) -htLep3t1minp = ROOT.TH1F("htlep3t1minp", "H_{T}^{Lep}", 300, 0., 600.) -htLep3t1topm = ROOT.TH1F("htlep3t1topm", "H_{T}^{Lep}", 300, 0., 600.) -htLep3t1topmhighptlep = ROOT.TH1F("htlep3t1topmhighptlep", "H_{T}^{Lep}", 300, 0., 600.) -htLep3w = ROOT.TH1F("htlep3w", "H_{T}^{Lep}", 300, 0., 600.) -htLep3top = ROOT.TH1F("htlep3top", "H_{T}^{Lep}", 300, 0., 600.) -htLep0 = ROOT.TH1F("htLep0", "H_{T}^{Lep}", 300, 0., 600.) -htLep1 = ROOT.TH1F("htLep1", "H_{T}^{Lep}", 300, 0., 600.) -htLep2 = ROOT.TH1F("htLep2", "H_{T}^{Lep}", 300, 0., 600.) -htLep3 = ROOT.TH1F("htLep3", "H_{T}^{Lep}", 300, 0., 600.) -ptJet0 = ROOT.TH1F("ptJet0", "p_{T} Of Leading Jet", 300, 0., 600.) - -topTagMassHistpremass= ROOT.TH1F("topTagMassHistpremass", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topTagugMassHistpremass= ROOT.TH1F("topTagugMassHistpremass", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topTagMassHist= ROOT.TH1F("topTagMassHist", "Mass of Top Candidate from Hadronic Jets type 1;Mass;Number", 300, 0., 600. ) -topTagMassHistPostTau32= ROOT.TH1F("topTagMassHistPostTau32", "Mass of Top Candidate from Hadronic Jets type 1;Mass;Number", 300, 0., 600. ) -topTagMassHistPostBDMax= ROOT.TH1F("topTagMassHistPostBDMax", "Mass of Top Candidate from Hadronic Jets type 1;Mass;Number", 300, 0., 600. ) - -topTagptHist= ROOT.TH1F("topTagptHist", "Pt of Top Candidate from Hadronic Jets type 1;Mass;Number", 375, 0., 1500. ) -topTagtau32Hist= ROOT.TH1F("topTagtau32Hist", "tau32 of Top Candidate from Hadronic Jets type 1;Tau32;Number", 150, 0., 1.5 ) -topTagBMaxHist= ROOT.TH1F("topTagBMaxHist", "BMax of Top Candidate from Hadronic Jets type 1;CSV;Number", 150, 0., 1.5 ) - - -topTagugmtau32Hist= ROOT.TH1F("topTagugmtau32Hist", "tau32 of Top Candidate from Hadronic Jets type 1;Tau32;Number", 150, 0., 1.5 ) -topTagugmBMaxHist= ROOT.TH1F("topTagugmBMaxHist", "BMax of Top Candidate from Hadronic Jets type 1;CSV;Number", 150, 0., 1.5 ) - - -ptlep= ROOT.TH1F("ptlep", "Pt Leptonic top", 150, 0., 1500. ) -topTagptHistprecuts= ROOT.TH1F("topTagptHistprecuts", "Pt of Top Candidate from Hadronic Jets type 1;Mass;Number", 375, 0., 1500. ) -topTagptHistprept= ROOT.TH1F("topTagptHistprept", "Pt of Top Candidate from Hadronic Jets type 1;Mass;Number", 375, 0., 1500. ) -minPairHist= ROOT.TH1F("minPairHist", "Minimum Pairwise mass", 150, 0., 150. ) -WMassPairHist= ROOT.TH1F("WMassPairHist", "Minimum Pairwise mass", 150, 0., 150. ) -WMassPairHistPtrel= ROOT.TH1F("WMassPairHistPtrel", "Minimum Pairwise mass", 150, 0., 150. ) -WMassPairHistDeltaR= ROOT.TH1F("WMassPairHistDeltaR", "Minimum Pairwise mass", 150, 0., 150. ) -WMassPairHistmu= ROOT.TH1F("WMassPairHistmu", "Minimum Pairwise mass", 150, 0., 150. ) -type1muHist = ROOT.TH1F("type1muHist", "#mu", 150, 0, 1.0) - -nvtxvsnsj = ROOT.TH2F("nvtxvsnsj","lep top pt vs number of subjets",40,0,80,11, -0.5, 10.5 ) -nvtxvsminmass = ROOT.TH2F("nvtxvsminmass","lep top pt vs number of subjets",40,0,80,150, 0., 150. ) -nvtxvstopmass = ROOT.TH2F("nvtxvstopmass","lep top pt vs number of subjets",40,0,80,300, 0., 600. ) -nvtxvstau32 = ROOT.TH2F("nvtxvstau32","lep top pt vs number of subjets",40,0,80,150, 0., 1.5 ) -nvtxvsbmax = ROOT.TH2F("nvtxvsbmax","lep top pt vs number of subjets",40,0,80,150, 0., 1.5 ) - -nsj= ROOT.TH1F("nsj", "number of subjets", 11, -0.5, 10.5 ) -nsjhighpt= ROOT.TH1F("nsjhighpt", "number of subjets", 11, -0.5, 10.5 ) -nsjmidptlep= ROOT.TH1F("nsjmidptlep", "number of subjets ptLep mid", 11, -0.5, 10.5 ) -nsjhighptlep= ROOT.TH1F("nsjhighptlep", "number of subjets ptLep high", 11, -0.5, 10.5 ) - -topcandmassprekin = ROOT.TH1F("topcandmassprekin", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topcandmasspostkin = ROOT.TH1F("topcandmasspostkin", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topcandmasspostnsj = ROOT.TH1F("topcandmasspostnsj", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topcandmasspostminmass = ROOT.TH1F("topcandmasspostminmass", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topcandmassposttau32 = ROOT.TH1F("topcandmassposttau32", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topcandmasspostbmax = ROOT.TH1F("topcandmasspostbmax", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) - - -topcandugmassprekin = ROOT.TH1F("topcandugmassprekin", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topcandugmasspostkin = ROOT.TH1F("topcandugmasspostkin", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topcandugmasspostnsj = ROOT.TH1F("topcandugmasspostnsj", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topcandugmasspostminmass = ROOT.TH1F("topcandugmasspostminmass", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topcandugmassposttau32 = ROOT.TH1F("topcandugmassposttau32", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topcandugmasspostbmax = ROOT.TH1F("topcandugmasspostbmax", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) - - - - -topTagMassHistpremasshighpt= ROOT.TH1F("topTagMassHistpremasshighpt", "Mass of Top Candidate from Hadronic Jets type 1 pre mass cut;Mass;Number", 300, 0., 600. ) -topTagMassHisthighpt= ROOT.TH1F("topTagMassHisthighpt", "Mass of Top Candidate from Hadronic Jets type 1;Mass;Number", 300, 0., 600. ) -topTagptHisthighpt= ROOT.TH1F("topTagptHisthighpt", "Pt of Top Candidate from Hadronic Jets type 1;Mass;Number", 375, 0., 1500. ) -minPairHisthighpt= ROOT.TH1F("minPairHisthighpt", "Minimum Pairwise mass", 150, 0., 150. ) - -ht3 = ROOT.TH1F("ht3", "HT", 200, 0., 2000.) -ht4 = ROOT.TH1F("ht4", "HT", 200, 0., 2000.) - -muHist = ROOT.TH1F("muHist", "#mu", 300, 0, 1.0) -muHisthighpt = ROOT.TH1F("muHisthighpt", "#mu", 300, 0, 1.0) -dRWbHist = ROOT.TH1F("dRWbHist", "#Delta R (W,b)", 300, 0.0, 5.0) - -mWCand = ROOT.TH1F("mWCand", "Mass of W Candidate from Hadronic Jets;Mass;Number", 200, 0., 200. ) -mWCandhighpt = ROOT.TH1F("mWCandhighpt", "Mass of W Candidate from Hadronic Jets;Mass;Number", 200, 0., 200. ) -mTopCand = ROOT.TH1F("mTopCand", "Mass of Top Candidate from Hadronic Jets;Mass;Number", 300, 0., 600. ) -mTopCandhighpt = ROOT.TH1F("mTopCandhighpt", "Mass of Top Candidate from Hadronic Jets;Mass;Number", 300, 0., 600. ) - -leptopptvstopmassprekin = ROOT.TH2F("leptopptvstopmassprekin","lep top pt vs topmass",150,0,1500,300, 0., 600. ) -leptopptvstopmasspostkin = ROOT.TH2F("leptopptvstopmasspostkin","lep top pt vs topmass",150,0,1500,300, 0., 600. ) -leptopptvstopmasspostnsj = ROOT.TH2F("leptopptvstopmasspostnsj","lep top pt vs topmass",150,0,1500,300, 0., 600. ) -leptopptvstopmasspostminmass = ROOT.TH2F("leptopptvstopmasspostminmass","lep top pt vs topmass",150,0,1500,300, 0., 600. ) -leptopptvstopmassposttau32 = ROOT.TH2F("leptopptvstopmassposttau32","lep top pt vs topmass",150,0,1500,300, 0., 600. ) -leptopptvstopmasspostbmax = ROOT.TH2F("leptopptvstopmasspostbmax","lep top pt vs topmass",150,0,1500,300, 0., 600. ) - - - - -leptopptvsnsj = ROOT.TH2F("leptopptvsnsj","lep top pt vs number of subjets",150,500,2000,11, -0.5, 10.5 ) - -mWCandVsMuCut = ROOT.TH2F("mWCandVsMuCut", "Mass of W candidate versus #mu cut", 20, 0, 200, 10, 0, 1.0) -mWCandVsMTopCand = ROOT.TH2F("mWCandVsMTopCand","WCand+bJet Mass vs WCand mass",200,0.,200.,600,0.,600.) - -mWCandVsPtWCand = ROOT.TH2F("mWCandVsPtWCand", "Mass of W candidate versus p_{T} of W candidate", 100, 0, 1000, 20, 0, 200) -mWCandVsPtWCandMuCut = ROOT.TH2F("mWCandVsPtWCandMuCut", "Mass of W candidate versus p_{T} of W candidate", 100, 0, 1000, 20, 0, 200) - - - -if options.makeResponse == True : - response = ROOT.RooUnfoldResponse(50, 300., 1300., 50, 300., 1300.) - response.SetName('response_pt') - ptTopCand = ROOT.TH1F("ptTopCand", "WCand+bJet p_{T};p_{T} (GeV/c);Number", 250, 0., 1000.) - ptWFromTopCand = ROOT.TH1F("ptWFromTopCand", "WCand in Type-2 top cand p_{T};p_{T} (GeV/c);Number", 250, 0., 1000.) - ptbFromTopCand = ROOT.TH1F("ptbFromTopCand", "bCand in Type-2 top cand p_{T};p_{T} (GeV/c);Number", 250, 0., 1000.) - - -events = Events (files) - -postfix = "" -if options.useLoose : - postfix = "Loose" - - -puHandle = Handle("int") -puLabel = ( "pileup", "npvRealTrue" ) - -npvHandle = Handle("unsigned int") -npvLabel = ( "pileup", "npv" ) - - -nsubCA8Handle = Handle( "vector > >") -nsubCA8Label = ( "nsub", "CA8P4" ) - -TopTau2Handle = Handle( "std::vector" ) -TopTau2Label = ( "nsub" , "Tau2") - -TopTau3Handle = Handle( "std::vector" ) -TopTau3Label = ( "nsub" , "Tau3") - -ak5JetPtHandle = Handle( "std::vector" ) -ak5JetPtLabel = ( "pfShyftTupleJets" + postfix + "AK5", "pt" ) -ak5JetEtaHandle = Handle( "std::vector" ) -ak5JetEtaLabel = ( "pfShyftTupleJets" + postfix + "AK5", "eta" ) -ak5JetPhiHandle = Handle( "std::vector" ) -ak5JetPhiLabel = ( "pfShyftTupleJets" + postfix + "AK5", "phi" ) -ak5JetMassHandle = Handle( "std::vector" ) -ak5JetMassLabel = ( "pfShyftTupleJets" + postfix + "AK5", "mass" ) -ak5JetDa0MassHandle = Handle( "std::vector" ) -ak5JetDa0MassLabel = ( "pfShyftTupleJets" + postfix + "AK5", "da0Mass" ) -ak5JetDa1MassHandle = Handle( "std::vector" ) -ak5JetDa1MassLabel = ( "pfShyftTupleJets" + postfix + "AK5", "da1Mass" ) -ak5JetCSVHandle = Handle( "std::vector" ) -ak5JetCSVLabel = ( "pfShyftTupleJets" + postfix + "AK5", "csv" ) -ak5JetSecvtxMassHandle = Handle( "std::vector" ) -ak5JetSecvtxMassLabel = ( "pfShyftTupleJets" + postfix + "AK5", "secvtxMass" ) - -ca8PrunedJetPtHandle = Handle( "std::vector" ) -ca8PrunedJetPtLabel = ( "pfShyftTupleJets" + postfix + "CA8Pruned", "pt" ) -ca8PrunedJetEtaHandle = Handle( "std::vector" ) -ca8PrunedJetEtaLabel = ( "pfShyftTupleJets" + postfix + "CA8Pruned", "eta" ) -ca8PrunedJetPhiHandle = Handle( "std::vector" ) -ca8PrunedJetPhiLabel = ( "pfShyftTupleJets" + postfix + "CA8Pruned", "phi" ) -ca8PrunedJetMassHandle = Handle( "std::vector" ) -ca8PrunedJetMassLabel = ( "pfShyftTupleJets" + postfix + "CA8Pruned", "mass" ) -ca8PrunedJetDa0MassHandle = Handle( "std::vector" ) -ca8PrunedJetDa0MassLabel = ( "pfShyftTupleJets" + postfix + "CA8Pruned", "da0Mass" ) -ca8PrunedJetDa1MassHandle = Handle( "std::vector" ) -ca8PrunedJetDa1MassLabel = ( "pfShyftTupleJets" + postfix + "CA8Pruned", "da1Mass" ) -ca8PrunedJetCSVHandle = Handle( "std::vector" ) -ca8PrunedJetCSVLabel = ( "pfShyftTupleJetsCA8" + postfix, "csv" ) - -topTagPtHandle = Handle( "std::vector" ) -topTagPtLabel = ( "pfShyftTupleJetsLooseTopTag" , "pt" ) -topTagEtaHandle = Handle( "std::vector" ) -topTagEtaLabel = ( "pfShyftTupleJetsLooseTopTag" , "eta" ) -topTagPhiHandle = Handle( "std::vector" ) -topTagPhiLabel = ( "pfShyftTupleJetsLooseTopTag" , "phi" ) -topTagMassHandle = Handle( "std::vector" ) -topTagMassLabel = ( "pfShyftTupleJetsLooseTopTag" , "mass" ) -topTagMinMassHandle = Handle( "std::vector" ) -topTagMinMassLabel = ( "pfShyftTupleJetsLooseTopTag" , "minMass" ) -topTagNSubjetsHandle = Handle( "std::vector" ) -topTagNSubjetsLabel = ( "pfShyftTupleJetsLooseTopTag" , "nSubjets" ) - -topTagsj0csvLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj0csv" ) -topTagsj0csvHandle = Handle( "std::vector" ) -topTagsj1csvLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj1csv" ) -topTagsj1csvHandle = Handle( "std::vector" ) -topTagsj2csvLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj2csv" ) -topTagsj2csvHandle = Handle( "std::vector" ) -topTagsj3csvLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj3csv" ) -topTagsj3csvHandle = Handle( "std::vector" ) - -topTagsj0ptLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj0pt" ) -topTagsj0ptHandle = Handle( "std::vector" ) -topTagsj1ptLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj1pt" ) -topTagsj1ptHandle = Handle( "std::vector" ) -topTagsj2ptLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj2pt" ) -topTagsj2ptHandle = Handle( "std::vector" ) -topTagsj3ptLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj3pt" ) -topTagsj3ptHandle = Handle( "std::vector" ) - -topTagsj0etaLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj0eta" ) -topTagsj0etaHandle = Handle( "std::vector" ) -topTagsj1etaLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj1eta" ) -topTagsj1etaHandle = Handle( "std::vector" ) -topTagsj2etaLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj2eta" ) -topTagsj2etaHandle = Handle( "std::vector" ) -topTagsj3etaLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj3eta" ) -topTagsj3etaHandle = Handle( "std::vector" ) - -topTagsj0phiLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj0phi" ) -topTagsj0phiHandle = Handle( "std::vector" ) -topTagsj1phiLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj1phi" ) -topTagsj1phiHandle = Handle( "std::vector" ) -topTagsj2phiLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj2phi" ) -topTagsj2phiHandle = Handle( "std::vector" ) -topTagsj3phiLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj3phi" ) -topTagsj3phiHandle = Handle( "std::vector" ) - -topTagsj0massLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj0mass" ) -topTagsj0massHandle = Handle( "std::vector" ) -topTagsj1massLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj1mass" ) -topTagsj1massHandle = Handle( "std::vector" ) -topTagsj2massLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj2mass" ) -topTagsj2massHandle = Handle( "std::vector" ) -topTagsj3massLabel = ( "pfShyftTupleJetsLooseTopTag" , "topsj3mass" ) -topTagsj3massHandle = Handle( "std::vector" ) - -muonPtHandle = Handle( "std::vector" ) -muonPtLabel = ( "pfShyftTupleMuons" + postfix, "pt" ) -muonEtaHandle = Handle( "std::vector" ) -muonEtaLabel = ( "pfShyftTupleMuons" + postfix, "eta" ) -muonPhiHandle = Handle( "std::vector" ) -muonPhiLabel = ( "pfShyftTupleMuons" + postfix, "phi" ) -muonPfisoHandle = Handle( "std::vector" ) -muonPfisoLabel = ( "pfShyftTupleMuons" + postfix, "pfisoPU" ) - -electronPtHandle = Handle( "std::vector" ) -electronPtLabel = ( "pfShyftTupleElectrons" + postfix, "pt" ) -electronEtaHandle = Handle( "std::vector" ) -electronEtaLabel = ( "pfShyftTupleElectrons" + postfix, "eta" ) -electronPhiHandle = Handle( "std::vector" ) -electronPhiLabel = ( "pfShyftTupleElectrons" + postfix, "phi" ) -electronPfisoHandle = Handle( "std::vector" ) -electronPfisoLabel = ( "pfShyftTupleElectrons" + postfix, "pfisoPU" ) - -metHandle = Handle( "std::vector" ) -metLabel = ("pfShyftTupleMET" + postfix, "pt" ) - -metphiHandle = Handle( "std::vector" ) -metphiLabel = ("pfShyftTupleMET" + postfix, "phi" ) - -if options.makeResponse == True : - genParticlesPtHandle = Handle( "std::vector") - genParticlesPtLabel = ( "pfShyftTupleGenParticles", "pt") - genParticlesEtaHandle = Handle( "std::vector") - genParticlesEtaLabel = ( "pfShyftTupleGenParticles", "eta") - genParticlesPhiHandle = Handle( "std::vector") - genParticlesPhiLabel = ( "pfShyftTupleGenParticles", "phi") - genParticlesMassHandle = Handle( "std::vector") - genParticlesMassLabel = ( "pfShyftTupleGenParticles", "mass") - genParticlesPdgIdHandle = Handle( "std::vector") - genParticlesPdgIdLabel = ( "pfShyftTupleGenParticles", "pdgId") - genParticlesStatusHandle = Handle( "std::vector") - genParticlesStatusLabel = ( "pfShyftTupleGenParticles", "status") - -goodEvents = [] -goodEventst1 = [] -mptv=0 -# loop over events -count = 0 -print "Start looping" -nhptbj=0 - -if options.set == 'mcatnlo': - print "Using MC@NLO Weights" - weightFile = ROOT.TFile("ptlepNewSelmcatnlo_weight.root") -elif options.set == 'powheg': - print "Using Powheg Weights" - weightFile = ROOT.TFile("ptlepNewSelpowheg_weight.root") -else: - print "Using MadGraph Weights" - weightFile = ROOT.TFile("ptlepNewSel_weight.root") -weightPlot = weightFile.Get("lepptweight") -if options.ptWeight == False : - print "Turning pt reweighting off" -for event in events: - weight = 1.0 - count = count + 1 - if count % 10000 == 0 : - print '--------- Processing Event ' + str(count) - if options.num != 'all': - if not (eventsbegin[ifile-1] <= count <= eventsend[ifile-1]): - continue - # if count > 900000: - #break - # if mptv % 100 == 0 : - # print '--------- mptv fail ' + str(mptv) - - - - # Find the top and antitop quarks. - # We also need to find the decay mode of the top and antitop quarks. - # To do so, we look for leptons, and use their charge to assign - # the correct decay mode to the correct quark. - topQuarks = [] - hadTop = None - lepTop = None - isSemiLeptonicGen = True - if options.makeResponse == True : - event.getByLabel( genParticlesPtLabel, genParticlesPtHandle ) - event.getByLabel( genParticlesEtaLabel, genParticlesEtaHandle ) - event.getByLabel( genParticlesPhiLabel, genParticlesPhiHandle ) - event.getByLabel( genParticlesMassLabel, genParticlesMassHandle ) - event.getByLabel( genParticlesPdgIdLabel, genParticlesPdgIdHandle ) - event.getByLabel( genParticlesStatusLabel, genParticlesStatusHandle ) - - genParticlesPt = genParticlesPtHandle.product() - genParticlesEta = genParticlesEtaHandle.product() - genParticlesPhi = genParticlesPhiHandle.product() - genParticlesMass = genParticlesMassHandle.product() - genParticlesPdgId = genParticlesPdgIdHandle.product() - genParticlesStatus = genParticlesStatusHandle.product() - - #print '------------' - p4Top = ROOT.TLorentzVector() - p4Antitop = ROOT.TLorentzVector() - topDecay = 0 # 0 = hadronic, 1 = leptonic - antitopDecay = 0 # 0 = hadronic, 1 = leptonic - - - - for igen in xrange( len(genParticlesPt) ) : - #print '{0:6.0f} {1:6.0f} {2:6.2f}'.format( genParticlesPdgId[igen], - # genParticlesStatus[igen], - # genParticlesPt[igen] ) - # Figure out the top vs antitop from charge of decay products. - # Need to know which is the hadronic one (top or antitop) to do unfolding - - if genParticlesStatus[igen] != 3 : - continue - if abs(genParticlesPdgId[igen]) < 6 : - continue - if abs(genParticlesPdgId[igen]) > 16 : - continue - - - if genParticlesStatus[igen] == 3 and genParticlesPdgId[igen] == 6 : - gen = ROOT.TLorentzVector() - gen.SetPtEtaPhiM( genParticlesPt[igen], - genParticlesEta[igen], - genParticlesPhi[igen], - genParticlesMass[igen] - ) - p4Top = gen - elif genParticlesStatus[igen] == 3 and genParticlesPdgId[igen] == -6 : - gen = ROOT.TLorentzVector() - gen.SetPtEtaPhiM( genParticlesPt[igen], - genParticlesEta[igen], - genParticlesPhi[igen], - genParticlesMass[igen] - ) - p4Antitop = gen - # If there is a lepton (e-, mu-, tau-) then the top is leptonic - elif genParticlesStatus[igen] == 3 and \ - ( genParticlesPdgId[igen] == -11 or genParticlesPdgId[igen] == -13 or genParticlesPdgId[igen] == -15) : - topDecay = 1 - # If there is an antilepton (e+, mu+, tau+) then the antitop is leptonic - elif genParticlesStatus[igen] == 3 and \ - ( genParticlesPdgId[igen] == 11 or genParticlesPdgId[igen] == 13 or genParticlesPdgId[igen] == 15) : - antitopDecay = 1 - - #print 'I think the top quark decay mode is ' + str(topDecay) - topQuarks.append( GenTopQuark( 6, p4Top, topDecay) ) - #print 'I think the top antiquark decay mode is ' + str(antitopDecay) - topQuarks.append( GenTopQuark( -6, p4Antitop, antitopDecay) ) - if topDecay + antitopDecay == 1 : - isSemiLeptonicGen = True - else : - isSemiLeptonicGen = False - - # If we are filling the response matrix, don't - # consider "volunteer" events that pass the selection - # even though they aren't really semileptonic events. - if isSemiLeptonicGen == False : - continue - - if topDecay == 0 : - hadTop = topQuarks[0] - lepTop = topQuarks[1] - else : - hadTop = topQuarks[1] - lepTop = topQuarks[0] - - - lepType = 0 # Let 0 = muon, 1 = electron - event.getByLabel (muonPtLabel, muonPtHandle) - if not muonPtHandle.isValid(): - mptv+=1 - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - muonPts = muonPtHandle.product() - - if True : - event.getByLabel (muonPfisoLabel, muonPfisoHandle) - if not muonPfisoHandle.isValid(): - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - muonPfisos = muonPfisoHandle.product() - - nMuonsVal = 0 - for imuonPt in range(0,len(muonPts)): - muonPt = muonPts[imuonPt] - if muonPt > 45.0 : - if options.useLoose : -# print 'imu = ' + str(imuonPt) + ', iso/pt = ' + str( muonPfisos[imuonPt] ) + '/' + str(muonPt) + ' = ' + str(muonPfisos[imuonPt]/muonPt) - if muonPfisos[imuonPt] / muonPt < 0.12 : - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - else : - nMuonsVal += 1 - lepType = 0 - else : - if muonPfisos[imuonPt] / muonPt > 0.12 : - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - nMuonsVal += 1 - lepType = 0 - if options.pileup != 'none': - event.getByLabel (puLabel, puHandle) - PileUp = puHandle.product() - bin1 = PilePlot.FindBin(PileUp[0]) - weight *= PilePlot.GetBinContent(bin1) - #print "weight is " + str(weight) - - event.getByLabel (npvLabel, npvHandle) - numvert = npvHandle.product() - - - nMuons.Fill( nMuonsVal,weight ) - #print "Filling???" - if nMuonsVal > 0 : - ptMu.Fill( muonPts[0],weight ) - - event.getByLabel (electronPtLabel, electronPtHandle) - electronPts = electronPtHandle.product() - - event.getByLabel (topTagEtaLabel, topTagEtaHandle) - topTagEta = topTagEtaHandle.product() - if not options.type2: - if len(topTagEta)<1: - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - if options.etabin != 'all': - if options.etabin == 'range1': - if abs(topTagEta[0])>1.0: - continue - if options.etabin == 'range2': - if abs(topTagEta[0])<=1.0: - continue - - nElectronsVal = 0 - # if nMuonsVal == 0 : - # for ielectronPt in range(0,len(electronPts)): - # electronPt = electronPts[ielectronPt] - # if electronPt > 45.0 : - # if options.useLoose : - # if electronPfisos[ielectronPt] / electronPt < 0.2 : - # continue - # else : - # nElectronsVal += 1 - # lepType = 1 - # else : - # nElectronsVal += 1 - # lepType = 1 - # - nElectrons.Fill( nElectronsVal,weight ) - - - if not options.muOnly : - if nMuonsVal + nElectronsVal != 1 : - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - else : - if nMuonsVal != 1: - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - - - # Now look at the rest of the lepton information. - # We will classify jets based on hemispheres, defined - # by the lepton. - - lepMass = 0.0 - if lepType == 0 : - event.getByLabel (muonEtaLabel, muonEtaHandle) - muonEtas = muonEtaHandle.product() - event.getByLabel (muonPhiLabel, muonPhiHandle) - muonPhis = muonPhiHandle.product() - - lepPt = muonPts[0] - lepEta = muonEtas[0] - lepPhi = muonPhis[0] - lepMass = 0.105 - else : - event.getByLabel (electronEtaLabel, electronEtaHandle) - electronEtas = electronEtaHandle.product() - event.getByLabel (electronPhiLabel, electronPhiHandle) - electronPhis = electronPhiHandle.product() - - lepPt = electronPts[0] - lepEta = electronEtas[0] - lepPhi = electronPhis[0] - - lepP4 = ROOT.TLorentzVector() - lepP4.SetPtEtaPhiM( lepPt, lepEta, lepPhi, lepMass ) - - - event.getByLabel (metLabel, metHandle) - mets = metHandle.product() - met = mets[0] - - htLepVal = met + lepP4.Perp() - - htLep0.Fill( htLepVal,weight ) - ptMET0.Fill( met,weight ) - - event.getByLabel (ak5JetPtLabel, ak5JetPtHandle) - ak5JetPts = ak5JetPtHandle.product() - - - nJetsVal = 0 - for jetPt in ak5JetPts: - if jetPt > 30.0 : - nJetsVal += 1 - - nJets.Fill( nJetsVal,weight ) - - # Require >= 2 AK5 jets above 30 GeV - if nJetsVal < 2 : - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - - htLep1.Fill( htLepVal,weight ) - ptMET1.Fill( met,weight ) - ptJet0.Fill( ak5JetPts[0],weight ) - - # Require leading jet pt to be pt > 200 GeV - - - htLep2.Fill( htLepVal,weight ) - ptMET2.Fill( met,weight ) - event.getByLabel (ak5JetCSVLabel, ak5JetCSVHandle) - ak5JetCSVs = ak5JetCSVHandle.product() - - event.getByLabel (ak5JetSecvtxMassLabel, ak5JetSecvtxMassHandle) - ak5JetSecvtxMasses = ak5JetSecvtxMassHandle.product() - - - ntags = 0 - for csv in ak5JetCSVs : - if csv > options.bDiscCut : - ntags += 1 - - if ntags < 1 : - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - - - - - htLep3.Fill( htLepVal,weight ) - ptMET3.Fill( met,weight ) - # Break the jets up by hemisphere - event.getByLabel (ak5JetEtaLabel, ak5JetEtaHandle) - ak5JetEtas = ak5JetEtaHandle.product() - event.getByLabel (ak5JetPhiLabel, ak5JetPhiHandle) - ak5JetPhis = ak5JetPhiHandle.product() - event.getByLabel (ak5JetMassLabel, ak5JetMassHandle) - ak5JetMasss = ak5JetMassHandle.product() - - ## event.getByLabel (ca8JetEtaLabel, ca8JetEtaHandle) - ## ca8JetEtas = ca8JetEtaHandle.product() - ## event.getByLabel (ca8JetPhiLabel, ca8JetPhiHandle) - ## ca8JetPhis = ca8JetPhiHandle.product() - ## event.getByLabel (ca8JetMassLabel, ca8JetMassHandle) - ## ca8JetMasss = ca8JetMassHandle.product() - ## event.getByLabel (ca8JetDa0MassLabel, ca8JetDa0MassHandle) - ## ca8JetDa0Masses = ca8JetDa0MassHandle.product() - ## event.getByLabel (ca8JetDa1MassLabel, ca8JetDa1MassHandle) - ## ca8JetDa1Masses = ca8JetDa1MassHandle.product() - - - event.getByLabel (topTagPtLabel, topTagPtHandle) - topTagPt = topTagPtHandle.product() - - event.getByLabel (topTagPhiLabel, topTagPhiHandle) - topTagPhi = topTagPhiHandle.product() - event.getByLabel (topTagEtaLabel, topTagEtaHandle) - topTagEta = topTagEtaHandle.product() - event.getByLabel (topTagMassLabel, topTagMassHandle) - topTagMass = topTagMassHandle.product() - event.getByLabel (topTagMinMassLabel, topTagMinMassHandle) - topTagMinMass = topTagMinMassHandle.product() - event.getByLabel (topTagNSubjetsLabel, topTagNSubjetsHandle) - topTagNSub = topTagNSubjetsHandle.product() - - - event.getByLabel (topTagsj0ptLabel, topTagsj0ptHandle) - Topsj0pt = topTagsj0ptHandle.product() - - event.getByLabel (topTagsj0etaLabel, topTagsj0etaHandle) - Topsj0eta = topTagsj0etaHandle.product() - - event.getByLabel (topTagsj0phiLabel, topTagsj0phiHandle) - Topsj0phi = topTagsj0phiHandle.product() - - event.getByLabel (topTagsj0massLabel, topTagsj0massHandle) - Topsj0mass = topTagsj0massHandle.product() - - - hadJets = [] - hadJetsMu = [] - hadJetsBDisc = [] - lepJets = [] - lepcsvs = [] - lepVtxMass = [] - ht = htLepVal - event.getByLabel (metphiLabel, metphiHandle) - metphis = metphiHandle.product() - metphi = metphis[0] - #print metphi - metv = ROOT.TLorentzVector() - metv.SetPtEtaPhiM( met, 0.0, metphi, 0.0) - - # Use AK5 jets for the leptonic side, and - # use CA8 jets for the hadronic side (below). This is - # because we want to fit the secondary vertex - # mass of the lepton-side b-jet. - for ijet in range(0,len(ak5JetPts)) : - if ak5JetPts[ijet] > 30.0 : - jet = ROOT.TLorentzVector() - jet.SetPtEtaPhiM( ak5JetPts[ijet], ak5JetEtas[ijet], ak5JetPhis[ijet], ak5JetMasss[ijet] ) - ht += jet.Perp() - if jet.DeltaR( lepP4 ) < ROOT.TMath.Pi() / 2.0 : - lepJets.append( jet ) - lepcsvs.append(ak5JetCSVs[ijet]) - lepVtxMass.append( ak5JetSecvtxMasses[ijet] ) - - - ht3.Fill( ht,weight ) - - - if options.htCut is not None and ht < options.htCut : - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - - - ## ht4.Fill( ht,weight ) - ## highestMassJetIndex = -1 - ## highestJetMass = -1.0 - ## bJetCandIndex = -1 - ## # Find highest mass jet for W-candidate - ## for ijet in range(0,len(hadJets)): - ## antitagged = hadJetsBDisc[ijet] < options.bDiscCut or options.bDiscCut < 0.0 - ## if hadJets[ijet].M() > highestJetMass and antitagged : - ## highestJetMass = hadJets[ijet].M() - ## highestMassJetIndex = ijet - - - - if len(lepJets) < 1: - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - - leptoppt = (metv+lepJets[0]+lepP4).Perp() - pt1 = leptoppt - Bin = weightPlot.FindBin(pt1) - ttweight = weightPlot.GetBinContent(Bin) - - - ntagslep=0 - - ugmass = -1 - - sjmass = [] - sjeta = [] - sjphi = [] - sjpt = [] - - - sjmass.append(Topsj0mass[0]) - sjpt.append(Topsj0pt[0]) - sjeta.append(Topsj0eta[0]) - sjphi.append(Topsj0phi[0]) - if topTagNSub[0] > 1: - event.getByLabel (topTagsj1ptLabel, topTagsj1ptHandle) - Topsj1pt = topTagsj1ptHandle.product() - - event.getByLabel (topTagsj1etaLabel, topTagsj1etaHandle) - Topsj1eta = topTagsj1etaHandle.product() - - event.getByLabel (topTagsj1phiLabel, topTagsj1phiHandle) - Topsj1phi = topTagsj1phiHandle.product() - - event.getByLabel (topTagsj1massLabel, topTagsj1massHandle) - Topsj1mass = topTagsj1massHandle.product() - - - sjmass.append(Topsj1mass[0]) - sjpt.append(Topsj1pt[0]) - sjeta.append(Topsj1eta[0]) - sjphi.append(Topsj1phi[0]) - if topTagNSub[0] > 2: - event.getByLabel (topTagsj2ptLabel, topTagsj2ptHandle) - Topsj2pt = topTagsj2ptHandle.product() - - - event.getByLabel (topTagsj2etaLabel, topTagsj2etaHandle) - Topsj2eta = topTagsj2etaHandle.product() - - event.getByLabel (topTagsj2phiLabel, topTagsj2phiHandle) - Topsj2phi = topTagsj2phiHandle.product() - - event.getByLabel (topTagsj2massLabel, topTagsj2massHandle) - Topsj2mass = topTagsj2massHandle.product() - - sjmass.append(Topsj2mass[0]) - sjpt.append(Topsj2pt[0]) - sjeta.append(Topsj2eta[0]) - sjphi.append(Topsj2phi[0]) - - if topTagNSub[0] > 3: - event.getByLabel (topTagsj3ptLabel, topTagsj3ptHandle) - Topsj3pt = topTagsj3ptHandle.product() - - event.getByLabel (topTagsj3etaLabel, topTagsj3etaHandle) - Topsj3eta = topTagsj3etaHandle.product() - - event.getByLabel (topTagsj3phiLabel, topTagsj3phiHandle) - Topsj3phi = topTagsj3phiHandle.product() - - event.getByLabel (topTagsj3massLabel, topTagsj3massHandle) - Topsj3mass = topTagsj3massHandle.product() - - sjmass.append(Topsj3mass[0]) - sjpt.append(Topsj3pt[0]) - sjeta.append(Topsj3eta[0]) - sjphi.append(Topsj3phi[0]) - - - - - sj = ROOT.TLorentzVector() - sjets = [] - sjets.append(ROOT.TLorentzVector()) - sjets[0].SetPtEtaPhiM( sjpt[0], sjeta[0], sjphi[0], sjmass[0] ) - topcomp = ROOT.TLorentzVector() - topcomp.SetPtEtaPhiM( sjpt[0], sjeta[0], sjphi[0], sjmass[0] ) - #print str(topTagNSub[0]) + " subjets" - for isub in range(1,int(topTagNSub[0])): - #print isub - sj.SetPtEtaPhiM( sjpt[isub], sjeta[isub], sjphi[isub], sjmass[isub] ) - sjets.append(ROOT.TLorentzVector()) - sjets[isub].SetPtEtaPhiM( sjpt[isub], sjeta[isub], sjphi[isub], sjmass[isub] ) - #print sj.M() - topcomp+=sj - ugmass = topcomp.M() - - - for lepcsv in lepcsvs : - if lepcsv > options.bDiscCut : - ntagslep += 1 - if ntagslep<1: - if options.makeResponse == True : - response.Miss( hadTop.p4.Perp(), weight ) - continue - - #print "weight " + str(weight) - nvtx =float(numvert[0]) - - if len(topTagPt) > 0: - if options.ptWeight == True : - t1weight=weight*ttweight - else: - t1weight=weight - #if leptoppt > 400.0: - jet1 = ROOT.TLorentzVector() - jet1.SetPtEtaPhiM( topTagPt[0], topTagEta[0], topTagPhi[0], topTagMass[0] ) - if ntagslep>0: - if (jet1.DeltaR( lepP4) > ROOT.TMath.Pi() / 2.0) : - ptlep.Fill(leptoppt,t1weight) - topcandmassprekin.Fill(topTagMass[0],t1weight) - - topcandugmassprekin.Fill(ugmass,t1weight) - leptopptvstopmassprekin.Fill(leptoppt,topTagMass[0],t1weight) - if (topTagPt[0] > 400.): - - topcandmasspostkin.Fill(topTagMass[0],t1weight) - topcandugmasspostkin.Fill(ugmass,t1weight) - leptopptvstopmasspostkin.Fill(leptoppt,topTagMass[0],t1weight) - if topTagNSub[0] > 2: - leptopptvstopmasspostnsj.Fill(leptoppt,topTagMass[0],t1weight) - topcandmasspostnsj.Fill(topTagMass[0],t1weight) - topcandugmasspostnsj.Fill(ugmass,t1weight) - if topTagMinMass[0] > 50. : - leptopptvstopmasspostminmass.Fill(leptoppt,topTagMass[0],t1weight) - topcandmasspostminmass.Fill(topTagMass[0],t1weight) - topcandugmasspostminmass.Fill(ugmass,t1weight) - event.getByLabel (nsubCA8Label, nsubCA8Handle) - nsubCA8Jets = nsubCA8Handle.product() - - event.getByLabel (topTagsj0csvLabel, topTagsj0csvHandle) - Topsj0BDiscCSV = topTagsj0csvHandle.product() - - event.getByLabel (topTagsj1csvLabel, topTagsj1csvHandle) - Topsj1BDiscCSV = topTagsj1csvHandle.product() - - event.getByLabel (topTagsj2csvLabel, topTagsj2csvHandle) - Topsj2BDiscCSV = topTagsj2csvHandle.product() - - - event.getByLabel (TopTau2Label, TopTau2Handle) - Tau2 = TopTau2Handle.product() - - event.getByLabel (TopTau3Label, TopTau3Handle) - Tau3 = TopTau3Handle.product() - - index = -1 - - for iCAjet in range(0,len(nsubCA8Jets)): - CAjetTLV = ROOT.TLorentzVector() - CAjetTLV.SetPtEtaPhiM( nsubCA8Jets[iCAjet].pt(), nsubCA8Jets[iCAjet].eta(), nsubCA8Jets[iCAjet].phi(), nsubCA8Jets[iCAjet].mass() ) - if (abs(jet1.DeltaR(CAjetTLV))<0.5): - index = iCAjet - break - - TauDisc = Tau3[index]/Tau2[index] - - if TauDisc<0.6: - leptopptvstopmassposttau32.Fill(leptoppt,topTagMass[0],t1weight) - topcandmassposttau32.Fill(topTagMass[0],t1weight) - topcandugmassposttau32.Fill(ugmass,t1weight) - BDMax = max(Topsj0BDiscCSV[0],Topsj1BDiscCSV[0],Topsj2BDiscCSV[0]) - if BDMax>0.679: - leptopptvstopmasspostbmax.Fill(leptoppt,topTagMass[0],t1weight) - topcandmasspostbmax.Fill(topTagMass[0],t1weight) - topcandugmasspostbmax.Fill(ugmass,t1weight) - - - - - #for ijet in range(0, len(topTagPt)) : - ijet=0 - jet = ROOT.TLorentzVector() - jet.SetPtEtaPhiM( topTagPt[ijet], topTagEta[ijet], topTagPhi[ijet], topTagMass[ijet] ) - if not len(lepJets) == 0: - if options.ptWeight == True : - t1weight=weight*ttweight - else: - t1weight=weight - - passSelection = False - if (jet.DeltaR( lepP4) > ROOT.TMath.Pi() / 2.0) : - if (topTagPt[ijet] > 250.): - topTagptHistprecuts.Fill(topTagPt[ijet],t1weight) - htLep3t1kin.Fill(htLepVal,t1weight) - nsj.Fill(topTagNSub[ijet],t1weight) - nvtxvsnsj.Fill(nvtx,topTagNSub[ijet],t1weight) - if topTagNSub[ijet] > 2: - minPairHist.Fill(topTagMinMass[ijet],t1weight) - nvtxvsminmass.Fill(nvtx,topTagMinMass[ijet],t1weight) - #WMassPairHist.Fill(WpairMass,t1weight) - #WMassPairHistPtrel.Fill(WpairMassptrel,t1weight) - #WMassPairHistDeltaR.Fill(WpairMassDeltaR,t1weight) - #WMassPairHistmu.Fill(WPairMassmu,t1weight) - #type1muHist.Fill(type1mu,t1weight) - if topTagMinMass[ijet] > 50. : - - htLep3t1minp.Fill(htLepVal,t1weight) - topTagMassHistpremass.Fill(topTagMass[ijet],t1weight) - nvtxvstopmass.Fill(nvtx,topTagMass[ijet],t1weight) - - topTagugMassHistpremass.Fill(topTagMass[ijet],t1weight) - - event.getByLabel (nsubCA8Label, nsubCA8Handle) - nsubCA8Jets = nsubCA8Handle.product() - - event.getByLabel (topTagsj0csvLabel, topTagsj0csvHandle) - Topsj0BDiscCSV = topTagsj0csvHandle.product() - - event.getByLabel (topTagsj1csvLabel, topTagsj1csvHandle) - Topsj1BDiscCSV = topTagsj1csvHandle.product() - - event.getByLabel (topTagsj2csvLabel, topTagsj2csvHandle) - Topsj2BDiscCSV = topTagsj2csvHandle.product() - - - event.getByLabel (TopTau2Label, TopTau2Handle) - Tau2 = TopTau2Handle.product() - - event.getByLabel (TopTau3Label, TopTau3Handle) - Tau3 = TopTau3Handle.product() - - index = -1 - - for iCAjet in range(0,len(nsubCA8Jets)): - CAjetTLV = ROOT.TLorentzVector() - CAjetTLV.SetPtEtaPhiM( nsubCA8Jets[iCAjet].pt(), nsubCA8Jets[iCAjet].eta(), nsubCA8Jets[iCAjet].phi(), nsubCA8Jets[iCAjet].mass() ) - if (abs(jet.DeltaR(CAjetTLV))<0.5): - index = iCAjet - break - - TauDisc = Tau3[index]/Tau2[index] - if ugmass > 150. and ugmass < 230.: - topTagugmtau32Hist.Fill(TauDisc,t1weight) - if TauDisc<0.6: - BDMax = max(Topsj0BDiscCSV[ijet],Topsj1BDiscCSV[ijet],Topsj2BDiscCSV[ijet]) - topTagugmBMaxHist.Fill(BDMax,t1weight) - - if topTagMass[ijet] > 140. and topTagMass[ijet] < 250.: - goodEventst1.append( [ event.object().id().run(), event.object().id().luminosityBlock(), event.object().id().event() ] ) - htLep3t1topm.Fill(htLepVal,t1weight) - topTagptHist.Fill(topTagPt[ijet],t1weight) - topTagMassHist.Fill(topTagMass[ijet],t1weight) - topTagtau32Hist.Fill(TauDisc,t1weight) - - nvtxvstau32.Fill(nvtx,TauDisc,t1weight) - - passSelection = True - - if TauDisc<0.6: - BDMax = max(Topsj0BDiscCSV[ijet],Topsj1BDiscCSV[ijet],Topsj2BDiscCSV[ijet]) - topTagBMaxHist.Fill(BDMax,t1weight) - nvtxvsbmax.Fill(nvtx,BDMax,t1weight) - topTagMassHistPostTau32.Fill(topTagMass[ijet],t1weight) - if BDMax>0.679: - topTagMassHistPostBDMax.Fill(topTagMass[ijet],t1weight) - if options.makeResponse == True : - if passSelection == True : - response.Fill( topTagPt[ijet], hadTop.p4.Perp(), t1weight ) - else : - response.Miss( hadTop.p4.Perp(), t1weight ) - - - -print 'Total Events: ' + str(count) -print 'isvalid() cuts: ' + str(mptv) -print 'percent cuts: ' + str((100*mptv)/count)+'%' -f.cd() - -if options.makeResponse : - response.Write() - -f.Write() - -f.Close() - -if options.printEvents : - Outf1 = open("type2skim/" + name + "Type2Skim.txt", "w") - sys.stdout = Outf1 - for goodEvent in goodEvents : - print '{0:12.0f}:{1:12.0f}:{2:12.0f}'.format( - goodEvent[0], goodEvent[1], goodEvent[2] - ) - Outf2 = open("type1skim/" + name + "Type1Skim.txt", "w") - sys.stdout = Outf2 - for goodEvent in goodEventst1 : - print '{0:12.0f}:{1:12.0f}:{2:12.0f}'.format( - goodEvent[0], goodEvent[1], goodEvent[2] - ) - - Outf3 = open("fullskim/" + name + "Skim.txt", "w") - sys.stdout = Outf3 - for goodEvent in goodEvents : - print '{0:12.0f}:{1:12.0f}:{2:12.0f}'.format( - goodEvent[0], goodEvent[1], goodEvent[2] - ) - for goodEvent in goodEventst1 : - print '{0:12.0f}:{1:12.0f}:{2:12.0f}'.format( - goodEvent[0], goodEvent[1], goodEvent[2] - ) - - -print "Total time = " + str( time.time() - start_time) + " seconds"