Merge pull request cms-sw#123 from MarkBaber/CMG_MiniAOD_Lite_V6_0_fr…

…om-CMSSW_7_0_6_AlphaT19Sep14 First commit of AlphaT configurations and code
gpetruc · Sep 22, 2014 · 4afffed · 4afffed
2 parents 50604c9 + be996ba
commit 4afffed
Show file tree

Hide file tree

Showing 8 changed files with 469 additions and 0 deletions.
diff --git a/CMGTools/RootTools/interface/AlphaT.h b/CMGTools/RootTools/interface/AlphaT.h
@@ -0,0 +1,26 @@
+#ifndef AlphaT_h
+#define AlphaT_h
+
+#include <cmath>
+#include <numeric>
+#include <vector>
+
+
+/*  AlphaT
+ *
+ *  Calculates the AlphaT event variable for a given input jet collection
+ *
+ *
+ *  08 Sept 2014 - Ported from ICF code (Mark Baber)
+ */
+
+
+struct AlphaT {
+
+  static double getAlphaT( const std::vector<double>& et,
+		    const std::vector<double>& px,
+		    const std::vector<double>& py );
+
+};
+
+#endif // AlphaT_h
diff --git a/CMGTools/RootTools/src/AlphaT.cc b/CMGTools/RootTools/src/AlphaT.cc
@@ -0,0 +1,34 @@
+#include "CMGTools/RootTools/interface/AlphaT.h"
+
+#include <cmath>
+#include <numeric>
+#include <vector>
+
+
+double AlphaT::getAlphaT( const std::vector<double>& et,
+			  const std::vector<double>& px,
+			  const std::vector<double>& py ){
+
+    // Momentum sums in transverse plane
+    const double sum_et = accumulate( et.begin(), et.end(), 0. );
+    const double sum_px = accumulate( px.begin(), px.end(), 0. );
+    const double sum_py = accumulate( py.begin(), py.end(), 0. );
+
+    // Minimum Delta Et for two pseudo-jets
+    double min_delta_sum_et = -1.;
+    for ( unsigned i=0; i < unsigned(1<<(et.size()-1)); i++ ) { //@@ iterate through different combinations
+      double delta_sum_et = 0.;
+      std::vector<bool> jet;
+      for ( unsigned j=0; j < et.size(); j++ ) { //@@ iterate through jets
+	delta_sum_et += et[j] * ( 1 - 2 * (int(i>>j)&1) ); 
+      }
+      if ( ( fabs(delta_sum_et) < min_delta_sum_et || min_delta_sum_et < 0. ) ) {
+	min_delta_sum_et = fabs(delta_sum_et);
+      }
+    }
+    if ( min_delta_sum_et < 0. ) { return 0.; }
+
+    // Alpha_T
+    return ( 0.5 * ( sum_et - min_delta_sum_et ) / sqrt( sum_et*sum_et - (sum_px*sum_px+sum_py*sum_py) ) );
+
+  }
diff --git a/CMGTools/RootTools/src/classes.h b/CMGTools/RootTools/src/classes.h
@@ -9,6 +9,7 @@
 #include "CMGTools/RootTools/interface/Davismt2.h"
 #include "CMGTools/RootTools/interface/mt2w_bisect.h"
 #include "CMGTools/RootTools/interface/Hemisphere.h"
+#include "CMGTools/RootTools/interface/AlphaT.h"
 
 namespace {
   namespace {
@@ -22,5 +23,6 @@ namespace {
     Hemisphere hemisphere(vector<float> px, vector<float> py, vector<float> pz, vector<float> E, int hemi_seed, int hemi_association);
     Davismt2 mt2;
     mt2w_bisect::mt2w mt2wlept;
+    AlphaT alphaT;
   }
 }
diff --git a/CMGTools/RootTools/src/classes_def.xml b/CMGTools/RootTools/src/classes_def.xml
@@ -10,4 +10,5 @@
   <class name="TriggerBitChecker" transient="true" />
   <class name="CMGMuonCleanerBySegmentsAlgo" transient="true" />
   <class name="EGammaMvaEleEstimatorFWLite" transient="true" />
+  <class name="AlphaT" transient="true" />
 </lcgdict>
diff --git a/CMGTools/TTHAnalysis/cfg/run_susyAlphaT_cfg.py b/CMGTools/TTHAnalysis/cfg/run_susyAlphaT_cfg.py
@@ -0,0 +1,209 @@
+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 *
+
+##------------------------------------------
+## Redefine analyzer parameters
+##------------------------------------------
+
+# Muons
+#------------------------------
+ttHLepAna.loose_muon_pt               = 10.,
+ttHLepAna.loose_muon_eta              = 2.5,
+ttHLepAna.loose_muon_dxy              = 0.5
+ttHLepAna.loose_muon_dz               = 1.0
+#ttHLepAna.loose_muon_relIso           = 0.15
+
+# Electrons
+#------------------------------
+ttHLepAna.loose_electron_id           = "POG_Cuts_ID_2012_Veto"
+ttHLepAna.loose_electron_pt           = 10
+ttHLepAna.loose_electron_eta          = 2.5
+ttHLepAna.loose_electron_dxy          = 0.5
+ttHLepAna.loose_electron_dz           = 0.
+# ttHLepAna.loose_electron_relIso       = 0.15
+# ttHLepAna.loose_electron_lostHits     = 999 # no cut
+# ttHLepAna.inclusive_electron_lostHits = 999 # no cut
+# ttHLepAna.ele_isoCorr                 = "deltaBeta"
+# ttHLepAna.ele_tightId                 = "Cuts_2012"
+
+# Photons
+#------------------------------
+ttHPhoAna.ptMin                        = 25,
+ttHPhoAna.epaMax                       = 2.5,
+
+# Taus 
+#------------------------------
+ttHTauAna.etaMax         = 2.3
+ttHTauAna.dxyMax         = 99999.
+ttHTauAna.dzMax          = 99999.
+ttHTauAna.vetoLeptons    = False
+ttHTauAna.vetoLeptonsPOG = True
+
+
+# Jets (for event variables do apply the jetID and not PUID yet)
+#------------------------------
+ttHJetAna.relaxJetId      = False
+ttHJetAna.doPuId          = False
+ttHJetAna.jetEta          = 5.
+ttHJetAna.jetEtaCentral   = 3.
+ttHJetAna.jetPt           = 50.
+ttHJetAna.recalibrateJets = False
+ttHJetAna.jetLepDR        = 0.4
+ttHJetMCAna.smearJets     = False
+
+
+# Energy sums
+#------------------------------
+# NOTE: Currently energy sums are calculated with 40 GeV jets (ttHCoreEventAnalyzer.py)
+#       However, the input collection is cleanjets which have a 50 GeV cut so this is a labeling problem
+
+ttHJetMETSkim.htCut       = ('htJet50j', 0)
+ttHJetMETSkim.mhtCut      = ('htJet40j', 0)
+ttHJetMETSkim.nBJet       = ('CSVM', 0, "jet.pt() > 50")     # require at least 0 jets passing CSVM and pt > 50
+
+##------------------------------------------
+##  ISOLATED TRACK
+##------------------------------------------
+
+# those are the cuts for the nonEMu
+ttHIsoTrackAna = cfg.Analyzer(
+            'ttHIsoTrackAnalyzer',
+#            candidates='cmgCandidates',
+#            candidatesTypes='std::vector<cmg::Candidate>',
+            candidates      ='packedPFCandidates',
+            candidatesTypes ='std::vector<pat::PackedCandidate>',
+            ptMin           = 5, ### for pion 
+            ptMinEMU        = 5, ### for EMU
+            dzMax           = 0.1,
+            #####
+            isoDR           = 0.3,
+            ptPartMin       = 0,
+            dzPartMax       = 0.1,
+            maxAbsIso       = 8,
+            #####
+            MaxIsoSum       = 0.1, ### unused
+            MaxIsoSumEMU    = 0.2, ### unused
+            doSecondVeto    = False
+            )
+
+
+##------------------------------------------
+##  ALPHAT VARIABLES
+##------------------------------------------
+
+# Tree Producer
+ttHAlphaTAna = cfg.Analyzer(
+            'ttHAlphaTVarAnalyzer'
+            )
+
+##------------------------------------------
+##  PRODUCER
+##------------------------------------------
+
+from CMGTools.TTHAnalysis.samples.samples_8TeV_v517 import triggers_RA1_Bulk, triggers_RA1_Prompt, triggers_RA1_Parked, triggers_RA1_Single_Mu, triggers_RA1_Photon, triggers_RA1_Muon
+
+# Tree Producer
+treeProducer = cfg.Analyzer(
+        'treeProducerSusyAlphaT',
+        vectorTree = True,
+        PDFWeights = PDFWeights,
+        triggerBits = {
+            'Bulk'     : triggers_RA1_Bulk,
+            'Prompt'   : triggers_RA1_Prompt,
+            'Parked'   : triggers_RA1_Parked,
+            'SingleMu' : triggers_RA1_Single_Mu,
+            'Photon'   : triggers_RA1_Photon,
+            'Muon'     : triggers_RA1_Muon,
+            }
+        )
+
+#-------- SAMPLES AND TRIGGERS -----------
+from CMGTools.TTHAnalysis.samples.samples_13TeV_CSA14 import *
+
+# CSA13 PU20bx25 samples: DYJetsM50_PU20bx25, DYJetsM50pythia6_PU20bx25, DYJetsM50_HT200to400_PU20bx25, DYJetsM50_HT400to600_PU20bx25, DYJetsM50_HT600toInf_PU20bx25, DYJetsMuMuM50_PtZ180_PU20bx25, DYJetsMuMuM6pythia8_PU20bx25, DYJetsMuMuM15pythia8_PU20bx25, DYJetsMuMuM50pythia8_PU20bx25, DYJetsEEpythia8_PU20bx25, DYJetsMuMupythia8_PU20bx25, EWKWmin_PU20bx25, EWKWplus_PU20bx25, EWKZjj_PU20bx25, EleGun_PU20bx25, GGHTauTau_PU20bx25, GGHZZ4L_PU20bx25, GJet_PU20bx25, JPsiPt20_PU20bx25, JPsiPt7_PU20bx25, MinBias_PU20bx25, MuMinGunPt100_PU20bx25, MuMinGunPt10_PU20bx25, MuPlusGunPt100_PU20bx25, MuPlusGunPt10_PU20bx25, NeutrinoGun_PU20bx25, QCDEM_20to30_PU20bx25, QCDEM_30to80_PU20bx25, QCDEM_80to170_PU20bx25, QCDMu_20to30_PU20bx25, QCDMu_30to50_PU20bx25, QCDMu_50to80_PU20bx25, QCDMu_80to120_PU20bx25, QCDMu_pythia6_120to170_PU20bx25, QCDMu_pythia6_20to30_PU20bx25, QCDMu_pythia6_30to50_PU20bx25, QCDMu_pythia6_50to80_PU20bx25, QCDMu_pythia6_80to120_PU20bx25, T1tttt_PU20bx25, TTHBB_PU20bx25, TTHGG_PU20bx25, TTHTauTau_PU20bx25, TTHWW_PU20bx25, TTHZZ4L_PU20bx25, TTJets_PU20bx25, TTJets_PUS14, TTpythia8_PU20bx25, VBFHBB_PU20bx25, VBFHGG_PU20bx25, VBFHWWSemi_PU20bx25, VBFHWW_PU20bx25, VBFHZG_PU20bx25, VBFHZZ4L_PU20bx25, VHMuMu_PU20bx25, VHTauTau_PU20bx25, VHWWInc_PU20bx25, VHWWLep_PU20bx25, VHZZ4L_PU20bx25, WENupyhia8_PU20bx25, WJets_PU20bx25, WminTau_PU20bx25, WplusMu_PU20bx25, WplusTau_PU20bx25, ZHBBInv_PU20bx25, ZHBBLL_PU20bx25, ZHLLInv_PU20bx25
+
+
+
+# Selected samples as defined on the AlphaT twiki
+WJetsToLNu   = [ WJetsToLNu_HT100to200_PU_S14_POSTLS170, WJetsToLNu_HT200to400_PU_S14_POSTLS170, WJetsToLNu_HT400to600_PU_S14_POSTLS170, WJetsToLNu_HT600toInf_PU_S14_POSTLS170]
+
+# Currently not defined in the samples file could be added from here: https://cmsweb.cern.ch/das/request?view=list&limit=100&instance=prod%2Fglobal&input=dataset%3D%2F*DYJetsToLL*13TeV*%2F*PU20bx25*%2F*AODSIM
+#DYJetsToLL  = []
+# Currently not defined in the samples file could be added from here: https://cmsweb.cern.ch/das/request?view=list&limit=100&instance=prod%2Fglobal&input=dataset%3D%2F*ZJetsToNuNu*13TeV*%2F*PU20bx25*%2F*AODSIM
+#ZJetsToNuNu = []
+# https://cmsweb.cern.ch/das/request?view=list&limit=100&instance=prod%2Fglobal&input=dataset%3D%2F*GJets*13TeV*%2F*PU20bx25*%2F*AODSIM
+#GJets       = []
+
+# NOT INCLUDING: /TTJets_MSDecaysCKM_central_Tune4C_13TeV-madgraph-tauola/Spring14miniaod-PU20bx25_POSTLS170_V5-v2/MINIAODSIM
+TTbar        = [ TTpythia8_PU20bx25 ]
+# https://cmsweb.cern.ch/das/request?view=list&limit=100&instance=prod%2Fglobal&input=dataset%3D%2FTToBLNu*13TeV*%2FSpring*PU20bx25*%2F*AODSIM
+#TToBLNu     = []
+# https://cmsweb.cern.ch/das/request?view=list&limit=100&instance=prod%2Fglobal&input=dataset%3D%2FSMS-T1qqqq*13TeV*%2FSpring*PU20bx25*%2F*AODSIM
+#T1qqqq       = []
+# https://cmsweb.cern.ch/das/request?view=list&limit=100&instance=prod%2Fglobal&input=dataset%3D%2FSMS-T1bbbb*13TeV*%2FSpring*PU20bx25*%2F*AODSIM
+#T1bbbb       = []
+T1tttt       = [ T1tttt_PU20bx25 ]
+
+
+#selectedComponents = [ SingleMu, DoubleElectron, TTHToWW_PUS14, DYJetsM50_PU20bx25, TTJets_PUS14 ]
+selectedComponents = []
+selectedComponents.extend( WJetsToLNu )
+selectedComponents.extend( TTbar )
+
+
+
+
+#-------- SEQUENCE
+
+sequence = cfg.Sequence(susyCoreSequence + [
+                        ttHIsoTrackAna,
+                        ttHAlphaTAna,
+                        treeProducer,
+                        ])
+
+
+#-------- HOW TO RUN
+test = 4
+
+# Test a single component, using a single thread.
+#--------------------------------------------------
+if test==1:
+    comp               = TTJets_PU20bx25
+    #comp.files = ['/afs/cern.ch/work/p/pandolf/CMSSW_7_0_6_patch1_2/src/CMGTools/TTHAnalysis/cfg/pickevents.root']
+    comp.files         = comp.files[:2]
+
+    selectedComponents = [comp]
+    comp.splitFactor   = 1
+#--------------------------------------------------
+
+# Test all components (1 thread per component).
+#--------------------------------------------------
+elif test==2:
+    for comp in selectedComponents:
+        comp.splitFactor = 1
+        comp.files       = comp.files[:1]
+#--------------------------------------------------
+
+# Run on local files
+#--------------------------------------------------
+elif test==4:
+    comp = TTJets_PU20bx25
+#    comp.name = 'TTJets'
+    #    comp.files = [ '/store/mc/Spring14miniaod/TT_Tune4C_13TeV-pythia8-tauola/MINIAODSIM/PU20bx25_POSTLS170_V5-v1/00000/063013AD-9907-E411-8135-0026189438BD.root' ]
+
+    comp.files = [ '/afs/cern.ch/user/m/mbaber/WORK/private/CSA14Samples/TT_Tune4C_13TeV-pythia8_PU20bx25.root' ]
+
+    selectedComponents = [comp]
+    comp.splitFactor = 1
+#--------------------------------------------------
+
+
+config = cfg.Config( components = selectedComponents,
+                     sequence = sequence )
+
+printComps(config.components, True)
+
diff --git a/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyAlphaT.py b/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyAlphaT.py
@@ -0,0 +1,93 @@
+from CMGTools.TTHAnalysis.analyzers.treeProducerSusyCore import *
+
+class treeProducerSusyAlphaT( treeProducerSusyCore ):
+
+    #-----------------------------------
+    # TREE PRODUCER FOR SUSY MULTILEPTONS
+    #-----------------------------------
+    def __init__(self, cfg_ana, cfg_comp, looperName):
+        super(treeProducerSusyAlphaT,self).__init__(cfg_ana, cfg_comp, looperName)
+
+        ## Declare what we want to fill (in addition to susy core ones)
+        self.globalVariables += [
+
+
+            #NTupleVariable("crossSection", lambda ev : ev.crossSection, help="process cross section in pb"),
+
+
+            # Energy sums
+            #----------------------------------------
+
+            NTupleVariable("ht",           lambda ev : ev.htJet40j,    help="H_{T} computed from only jets (with |eta|<2.5, pt > 40 GeV)"),
+            NTupleVariable("deltaPhiMin",  lambda ev : ev.deltaPhiMin, help="minimal deltaPhi between the MET and the four leading jets with pt>40 and eta<2.4"),
+            NTupleVariable("diffMetMht",   lambda ev : ev.diffMetMht,   help="abs( vec(mht) - vec(met) )"),
+            NTupleVariable("mht_pt",       lambda ev : ev.mhtJet40j,    help="H_{T}^{miss} computed from only jets (with |eta|<2.5, pt > 40 GeV)"),
+            NTupleVariable("mht_phi",      lambda ev : ev.mhtPhiJet40j, help="H_{T}^{miss} #phi computed from onlyy jets (with |eta|<2.5, pt > 40 GeV)"),
+            ##--------------------------------------------------
+
+            # Physics object multplicities
+            #----------------------------------------
+            NTupleVariable("nMuons10",     lambda ev: sum([l.pt() > 10 and abs(l.pdgId()) == 13 for l in ev.selectedLeptons]), int, help="Number of muons with pt > 10"),
+            NTupleVariable("nElectrons10", lambda ev: sum([l.pt() > 10 and abs(l.pdgId()) == 11 for l in ev.selectedLeptons]), int, help="Number of electrons with pt > 10"),
+            NTupleVariable("nTaus20",      lambda ev: sum([l.pt() > 20 for l in ev.selectedTaus]), int, help="Number of taus with pt > 20"),
+            NTupleVariable("nGammas20",    lambda ev: sum([l.pt() > 20 for l in ev.selectedPhotons]), int, help="Number of photons with pt > 20"),
+            NTupleVariable("nBJet40",      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"),
+
+            # AlphaT variables
+            #----------------------------------------
+
+            # NEED TO ADD:
+            # - M_{T}
+            # - metcorrected for photon, muon
+            # - mht/metcorrected
+
+            # NTupleVariable("nJet50",        lambda ev: len(ev.cleanJets), int, help="Number of jets with pt > 50, |eta|<3"),
+            # NTupleVariable("nJetFwd50",     lambda ev: len(ev.cleanJetsFwd), int, help="Number of jets with pt > 50, 3<|eta|<5"),
+            # NTupleVariable("nBJetMedium50", lambda ev: len(ev.bjetsMedium), int, help="Number of jets with pt > 50, |eta|<3 passing CSV medium"),
+            # NTupleVariable("htJet50",       lambda ev: ev.htJet50,     help="H_{T} computed from jets (with |eta|<3, pt > 50 GeV)"),
+            # NTupleVariable("mhtJet50",      lambda ev: ev.mhtJet50,    help="H_{T}^{miss} computed from jets (with |eta|<3, pt > 50 GeV)"),
+            # NTupleVariable("mht50OverMet",  lambda ev: ev.mht50OverMet,help="H_{T}^{miss}/E_{T}^{miss} computed from jets (with |eta|<3, pt > 50 GeV)"),
+
+
+            NTupleVariable("alphaT",        lambda ev: ev.alphaT, help="AlphaT computed using jets with pt > 50, |eta|<3"),
+
+            # # NTupleVariable("mT", lambda ev: ev.mT, help="M_{T} computed from lepton and type-1 corrected PFMET"),
+
+            # NTupleVariable("mZ1",           lambda ev: ev.bestZ1[0], help="Best m(ll) SF/OS"),
+            # NTupleVariable("pt2l",          lambda ev: ev.pt2l,      help="p_{T}(ll)"),
+            # NTupleVariable("m2l",           lambda ev: ev.m2l,       help="m(ll)"),
+
+            ##--------------------------------------------------
+            #            NTupleVariable("mtw", lambda ev: ev.mtw, int, help="mt(l,met)"),
+            #            NTupleVariable("mtwTau", lambda ev: ev.mtwTau, int, help="mt(tau,met)"),
+            #            NTupleVariable("IsoTrack_mtw", lambda ev: ev.mtwIsoTrack, int, help="mt(isoTrack,met)"),
+            ##--------------------------------------------------
+            #            NTupleVariable("minMWjj", lambda ev: ev.minMWjj, int, help="minMWjj"),
+            #            NTupleVariable("minMWjjPt", lambda ev: ev.minMWjjPt, int, help="minMWjjPt"),
+            #            NTupleVariable("bestMWjj", lambda ev: ev.bestMWjj, int, help="bestMWjj"),
+            #            NTupleVariable("bestMWjjPt", lambda ev: ev.bestMWjjPt, int, help="bestMWjjPt"),
+            #            NTupleVariable("bestMTopHad", lambda ev: ev.bestMTopHad, int, help="bestMTopHad"),
+            #            NTupleVariable("bestMTopHadPt", lambda ev: ev.bestMTopHadPt, int, help="bestMTopHadPt"),
+            ##--------------------------------------------------
+            ]
+
+        self.globalObjects.update({
+            # put more here
+            # "pseudoJet1"       : NTupleObject("pseudoJet1",     fourVectorType, help="pseudoJet1 for hemishphere"),
+            # "pseudoJet2"       : NTupleObject("pseudoJet2",     fourVectorType, help="pseudoJet2 for hemishphere"),
+            })
+        self.collections.update({
+            # put more here
+            "selectedLeptons"  : NTupleCollection("lep",      leptonTypeSusy,           50, help="Leptons after the preselection", filter=lambda l : l.pt()>10 ),
+            "selectedTaus"     : NTupleCollection("tau",      tauTypeSusy,              50, help="Taus after the preselection"),
+            "cleanJetsAll"     : NTupleCollection("jet",      jetTypeSusy,             100, help="all jets (w/ x-cleaning, w/ ID applied w/o PUID applied pt > 50 |eta| < 3) , sorted by pt", filter=lambda l : l.pt()>50  ),
+            "selectedPhotons"  : NTupleCollection("gamma",    photonTypeSusy,           50, help="photons with pt > 25 and loose cut based ID"),
+            "selectedIsoTrack" : NTupleCollection("isoTrack", isoTrackType,             50, help="isoTrack, sorted by pt"),
+            "genParticles"     : NTupleCollection("genPart",  genParticleWithMotherId, 200, help="all pruned genparticles"),
+            })
+
+        ## Book the variables, but only if we're called explicitly and not through a base class
+        if cfg_ana.name == "treeProducerSusyAlphaT":
+            self.initDone = True
+            self.declareVariables()
+