Skip to content

Instantly share code, notes, and snippets.

@tahuang1991
Created August 15, 2019 16:46
Show Gist options
  • Save tahuang1991/58d65307b573cb4ff58f115a08e301b9 to your computer and use it in GitHub Desktop.
Save tahuang1991/58d65307b573cb4ff58f115a08e301b9 to your computer and use it in GitHub Desktop.
## this file contains functions that check if a L1Mu is isolated
from Helpers import *
from hybridAlgorithmPtAssignment import *
import random
##_________________________________________________
def deltaPhi(phi1, phi2):
M_PI = 3.14159265358979323846264338328
dphi = phi1-phi2
if ( dphi > M_PI ):
dphi = dphi-M_PI*2.0
elif (dphi <= -M_PI):
dphi = dphi + M_PI*2.0
return dphi
def is_L1Mu_isolated(treeHits, L1Mu_index, vetoType):
"""checks if muon is isolated from track trigger track"""
L1Mu_L1Tk_dR_min = treeHits.L1Mu_L1Tk_dR_prop[L1Mu_index]
L1Mu_L1Tk_pt = treeHits.L1Mu_L1Tk_pt_prop[L1Mu_index]
verbose = False
if verbose:
print "Checking veto"
print "L1Mu_L1Tk_dR_min", L1Mu_L1Tk_dR_min
print "L1Mu_L1Tk_pt", L1Mu_L1Tk_pt
# The loose veto rejects prompt muons by matching a L1Tk within
# a radius R<0.12 with an L1Tk pT > 4 GeV. The medium and tight
# veto apply L1Tk pT cuts of 3 and 2 GeV respectively on L1Tk
# in R<0.12.
## loose
if vetoType == 1:
dR_largeCone = 0.12
ptCut_largeCone = 4
dR_smallCone = 0.12
ptCut_smallCone = 4
## medium
if vetoType == 2:
dR_largeCone = 0.12
ptCut_largeCone = 4
dR_smallCone = 0.12
ptCut_smallCone = 3
## tight
if vetoType == 3:
dR_largeCone = 0.12
ptCut_largeCone = 4
dR_smallCone = 0.12
ptCut_smallCone = 2
isMatched = False
isUnmatched = False
## check if matched or unmatched
## L1Tk should have a momentum above a certain threshold to be matched or unmatched
#if L1Mu_L1Tk_dR_min <= dR_largeCone and L1Mu_L1Tk_pt >= ptCut_largeCone: isUnmatched = True
if L1Mu_L1Tk_dR_min <= dR_smallCone and L1Mu_L1Tk_pt >= ptCut_smallCone: isMatched = True
## isolated means neither matched nor unmatched!
return (not isMatched) and (not isUnmatched)
##_________________________________________________
def isME11StubDisabled(failRate):
random_number = random.random()
return random_number < failRate
##_________________________________________________
def isME21StubDisabled(failRate):
random_number = random.random()
return random_number < failRate
##_________________________________________________
def fillDPhiHistogram( h_dphi_ME11_ME21, treeHits, etaCutMin, etaCutMax, npar, timesCharge = True, min_pt = 0, max_pt = 999,
stubCut=2,#not implemented yet
qualityCut=4,
hasME1Cut=False,
hasME2Cut=False):
## check if this event has L1Mus
if len(list(treeHits.L1Mu_pt))==0:
return max_prompt_L1Mu_pt, max_prompt_L1Mu_eta
pts = list(treeHits.L1Mu_pt)
#print "\tN L1Mus event", len(pts)
for i in range(0,len(pts)):
L1Mu_pt = treeHits.L1Mu_pt[i]
L1Mu_eta = treeHits.L1Mu_eta[i]
L1Mu_phi = treeHits.L1Mu_phi[i]
L1Mu_bx = treeHits.L1Mu_bx[i]
L1Mu_quality = treeHits.L1Mu_quality[i]
L1Mu_CSCTF_index = treeHits.L1Mu_CSCTF_index[i]
L1Mu_DTTF_index = treeHits.L1Mu_DTTF_index[i]
L1Mu_charge = treeHits.L1Mu_charge[i]
#print "L1Mu_charge", L1Mu_charge
if L1Mu_pt < min_pt or L1Mu_pt > max_pt: continue
## eta cut
if not (etaCutMin <= abs(L1Mu_eta) and abs(L1Mu_eta) <= etaCutMax): continue
## quality cut
if L1Mu_quality < qualityCut: continue
## BX cut
if abs(L1Mu_bx) != 0 and True: continue
## not a CSC muon
if L1Mu_CSCTF_index == -1: continue
has_CSC_ME1 = treeHits.CSCTF_phi1[L1Mu_CSCTF_index] != 99
has_CSC_ME2 = treeHits.CSCTF_phi2[L1Mu_CSCTF_index] != 99
if hasME1Cut and not has_CSC_ME1: continue
if hasME2Cut and not has_CSC_ME2: continue
## parity cut
CSCTF_ch1 = treeHits.CSCTF_ch1[L1Mu_CSCTF_index]
CSCTF_ch2 = treeHits.CSCTF_ch2[L1Mu_CSCTF_index]
CSCTF_isEven1 = CSCTF_ch1%2==0
CSCTF_isEven2 = CSCTF_ch2%2==0
if get_parity_ME11_ME21(CSCTF_isEven1, CSCTF_isEven2) != npar: continue
GE11_phi_L1 = treeHits.GE11_phi_L1[L1Mu_CSCTF_index]
GE11_phi_L2 = treeHits.GE11_phi_L2[L1Mu_CSCTF_index]
GE21_phi_L1 = treeHits.GE21_pad2_phi_L1[L1Mu_CSCTF_index]
GE21_phi_L2 = treeHits.GE21_pad2_phi_L2[L1Mu_CSCTF_index]
CSCTF_fit_phi1 = treeHits.CSCTF_fit_phi1[L1Mu_CSCTF_index]
CSCTF_fit_phi2 = treeHits.CSCTF_fit_phi2[L1Mu_CSCTF_index]
GE11_phi = getBestValue(GE11_phi_L1, GE11_phi_L2)
GE21_phi = getBestValue(GE21_phi_L1, GE21_phi_L2)
#GE11_ME11_dphi = - GE11_phi + CSCTF_fit_phi1
#GE21_ME21_dphi = - GE21_phi + CSCTF_fit_phi2
if CSCTF_fit_phi1 == 99 or CSCTF_fit_phi2 == 99:
continue
if GE11_phi == 99.:
GE11_ME11_dphi = 99#0.025
elif (timesCharge):
GE11_ME11_dphi = L1Mu_charge*deltaPhi(CSCTF_fit_phi1, GE11_phi)
else:
GE11_ME11_dphi = deltaPhi(CSCTF_fit_phi1, GE11_phi)
if GE21_phi == 99.:
GE21_ME21_dphi = 99#0.025
elif (timesCharge):
GE21_ME21_dphi = L1Mu_charge*deltaPhi(CSCTF_fit_phi2, GE21_phi)
else :
GE21_ME21_dphi = deltaPhi(CSCTF_fit_phi2, GE21_phi)
GE11_dphi_v2 = treeHits.CSCTF_gemdphi1[L1Mu_CSCTF_index]
GE21_dphi_v2 = treeHits.CSCTF_gemdphi2[L1Mu_CSCTF_index]
DPhi = abs(treeHits.CSCTF_L1_DPhi12_GE21[L1Mu_CSCTF_index])
"""
if ((npar==3 or npar ==1) and npar == get_parity_ME11_ME21(CSCTF_isEven1, CSCTF_isEven2)):
print "npar required ",npar, " get_parity ",get_parity_ME11_ME21(CSCTF_isEven1, CSCTF_isEven2)," GE11_phi ", GE11_phi," CSCTF_fit_phi1 ",CSCTF_fit_phi1," GE21_phi ", GE21_phi, " CSCTF_fit_phi2 ",CSCTF_fit_phi2," dphi_st2 ",GE21_ME21_dphi
print "another verioin treeHits.CSCTF_gemdphi1 ",GE11_dphi_v2," gemdphi2 ",GE21_dphi_v2," DPhi_dir ",DPhi,"\n"
print "GE11_phi_L1", GE11_phi_L1
print "GE11_phi_L2", GE11_phi_L2
print "GE21_phi_L1", GE21_phi_L1
print "GE21_phi_L2", GE21_phi_L2
print "ME11_phi", CSCTF_fit_phi1
print "ME21_phi", CSCTF_fit_phi2
print "GE11_phi", GE11_phi
print "GE21_phi", GE21_phi
print "GE11_ME11_dphi", GE11_ME11_dphi
print "GE21_ME21_dphi", GE21_ME21_dphi
print
"""
h_dphi_ME11_ME21.Fill(GE11_ME11_dphi, GE21_ME21_dphi)
##_________________________________________________
def fillPtEtaHistogram(ptHistogram,
etaHistogram,
treeHits,
doBXCut,
etaCutMin,
etaCutMax,
stubCut=2,
qualityCut=4,
hasME1Cut=False,
hasME2Cut=False,
hasME3Cut=False,
hasME4Cut=False,
hasME11Cut=False,
hasME21Cut=False,
hasGE11Cut=False,
hasGE21Cut=False,
hasME11ME21Cut=False,
hasGE11GE21Cut=False,
ME11FailRate=0,
ME21FailRate=0,
hasMB1Cut=False,
hasMB2Cut=False,
hasMB3Cut=False,
hasMB4Cut=False):
prompt_L1Mu_pt, prompt_L1Mu_eta, iL1Mu= getMaxPromptPtEtaL1MuEvent(treeHits,
doBXCut,
etaCutMin,
etaCutMax,
stubCut,
qualityCut,
hasME1Cut,
hasME2Cut,
hasME3Cut,
hasME4Cut,
hasME11Cut,
hasME21Cut,
hasGE11Cut,
hasGE21Cut,
hasME11ME21Cut,
hasGE11GE21Cut,
ME11FailRate,
ME21FailRate,
hasMB1Cut,
hasMB2Cut,
hasMB3Cut,
hasMB4Cut)
if (prompt_L1Mu_pt>0):
ptHistogram.Fill(prompt_L1Mu_pt)
## apply a 10 GeV pT cut for the eta histograms!!!
#if (prompt_L1Mu_pt>20):
#print "Max prompt pt, eta event", prompt_L1Mu_pt, prompt_L1Mu_eta
if (prompt_L1Mu_pt>=10):
etaHistogram.Fill(abs(prompt_L1Mu_eta))
def getMaxPromptPtEtaL1MuEvent(treeHits,
doBXCut,
etaCutMin,
etaCutMax,
stubCut,
qualityCut,
hasME1Cut=False,
hasME2Cut=False,
hasME3Cut=False,
hasME4Cut=False,
hasME11Cut=False,
hasME21Cut=False,
hasGE11Cut=False,
hasGE21Cut=False,
hasME11ME21Cut=False,
hasGE11GE21Cut=False,
ME11FailRate=0,
ME21FailRate=0,
hasMB1Cut=False,
hasMB2Cut=False,
hasMB3Cut=False,
hasMB4Cut=False):
max_prompt_L1Mu_pt = 0
max_prompt_L1Mu_eta = -99
iL1Mu = -1
## check if this event has L1Mus
if len(list(treeHits.L1Mu_pt))==0:
return max_prompt_L1Mu_pt, max_prompt_L1Mu_eta, iL1Mu
pts = list(treeHits.L1Mu_pt)
#print "\tN L1Mus event", len(pts)
for i in range(0,len(pts)):
L1Mu_pt = treeHits.L1Mu_pt[i]
L1Mu_eta = treeHits.L1Mu_eta[i]
L1Mu_phi = treeHits.L1Mu_phi[i]
L1Mu_bx = treeHits.L1Mu_bx[i]
L1Mu_quality = treeHits.L1Mu_quality[i]
L1Mu_CSCTF_index = treeHits.L1Mu_CSCTF_index[i]
L1Mu_DTTF_index = treeHits.L1Mu_DTTF_index[i]
L1Mu_eta_ME2 = -99
"""
## get the DTTF index
getDTTFindex = True
if getDTTFindex:
iDTTF_index = -1
iDTTF_dEta = 99
for iDTTF in range(0,len(treeHits.DTTF_bx)):
if treeHits.DTTF_bx[iDTTF] == L1Mu_bx and deltaPhi2(L1Mu_phi, treeHits.DTTF_phi[iDTTF])<0.001:
iDTTF_index = iDTTF
#print "found index",
L1Mu_DTTF_index = iDTTF_index
"""
## eta cut
if not (etaCutMin <= abs(L1Mu_eta) and abs(L1Mu_eta) <= etaCutMax): continue
## quality cut
if L1Mu_quality < qualityCut: continue
## BX cut
if abs(L1Mu_bx) != 0 and doBXCut: continue
## CSC quantities
has_CSC_ME1 = False
has_CSC_ME2 = False
has_actual_CSC_ME1 = False
has_actual_CSC_ME2 = False
has_CSC_ME3 = False
has_CSC_ME4 = False
has_CSC_ME11 = False
has_CSC_ME21 = False
is_CSC_ME11_disabled = False
is_CSC_ME21_disabled = False
has_GE11 = False
has_GE21 = False
GE11_dPhi = 99
CSC_ME1_ch = -1
CSC_ME2_ch = -2
nCSCStubs = 0
## DT quantities
has_DT_MB1 = False
has_DT_MB2 = False
has_DT_MB3 = False
has_DT_MB4 = False
nDTStubs = 0
#print L1Mu_CSCTF_index
#print L1Mu_CSCTF_index, len(treeHits.CSCTF_bx1)
if L1Mu_CSCTF_index != -1 and L1Mu_CSCTF_index < len(treeHits.CSCTF_phi1):
has_CSC_ME1 = treeHits.CSCTF_phi1[L1Mu_CSCTF_index] != 99
has_CSC_ME2 = treeHits.CSCTF_phi2[L1Mu_CSCTF_index] != 99
has_CSC_ME3 = treeHits.CSCTF_phi3[L1Mu_CSCTF_index] != 99
has_CSC_ME4 = treeHits.CSCTF_phi4[L1Mu_CSCTF_index] != 99
L1Mu_eta_ME2 = treeHits.CSCTF_eta2[L1Mu_CSCTF_index]
if hasME1Cut and not has_CSC_ME1: continue
if hasME2Cut and not has_CSC_ME2: continue
if hasME3Cut and not has_CSC_ME3: continue
if hasME4Cut and not has_CSC_ME4: continue
has_CSC_ME11 = treeHits.CSCTF_ri1[L1Mu_CSCTF_index] == 1 or treeHits.CSCTF_ri1[L1Mu_CSCTF_index] == 4
has_CSC_ME21 = treeHits.CSCTF_ri2[L1Mu_CSCTF_index] == 1
is_CSC_ME11_disabled = isME11StubDisabled(ME11FailRate)
is_CSC_ME21_disabled = isME21StubDisabled(ME21FailRate)
if hasME11Cut and not has_CSC_ME11 and not is_CSC_ME11_disabled: continue
if hasME21Cut and not has_CSC_ME21 and not is_CSC_ME21_disabled: continue
GE11_dPhi = treeHits.CSCTF_gemdphi1[L1Mu_CSCTF_index]
GE21_dPhi = treeHits.CSCTF_gemdphi2[L1Mu_CSCTF_index]
CSC_ME1_ch = treeHits.CSCTF_ch1[L1Mu_CSCTF_index]
CSC_ME2_ch = treeHits.CSCTF_ch2[L1Mu_CSCTF_index]
## check if stubs in station 1 and 2 can really be counted! -- they may be disabled
if not has_CSC_ME11: has_actual_CSC_ME1 = has_CSC_ME1
else: has_actual_CSC_ME1 = (not is_CSC_ME11_disabled)
if not has_CSC_ME21: has_actual_CSC_ME2 = has_CSC_ME2
else: has_actual_CSC_ME2 = (not is_CSC_ME21_disabled)
nCSCStubs = has_actual_CSC_ME1 + has_actual_CSC_ME2 + has_CSC_ME3 + has_CSC_ME4
if nCSCStubs < stubCut: continue
#print "is csc muon"
#print "\t", nCSCStubs
if (hasME11ME21Cut and not has_CSC_ME11 and not has_CSC_ME21): continue
## no bending angles
if is_CSC_ME11_disabled: GE11_dPhi = 99
if is_CSC_ME21_disabled: GE21_dPhi = 99
if hasGE11Cut and not passDPhicutTFTrack(1, CSC_ME1_ch, GE11_dPhi, L1Mu_pt):
#print "muon failed the GE11 cut", CSC_ME1_ch, GE11_dPhi, L1Mu_pt
continue
if hasGE21Cut and not passDPhicutTFTrack(2, CSC_ME2_ch, GE21_dPhi, L1Mu_pt):
#print "muon failed the GE21 cut", CSC_ME1_ch, GE11_dPhi, L1Mu_pt
continue
if hasGE11GE21Cut and GE11_dPhi == 99 and GE21_dPhi == 99 and iL1Mu >= 0: continue#only pass it if already got one L1Mu with bending
#if (hasGE11GE21Cut and (not passDPhicutTFTrack(1, CSC_ME1_ch, GE11_dPhi, L1Mu_pt)) and
# (not passDPhicutTFTrack(2, CSC_ME2_ch, GE21_dPhi, L1Mu_pt))): continue
if L1Mu_DTTF_index != -1 and L1Mu_DTTF_index < len(treeHits.DTTF_phi1):
has_DT_MB1 = treeHits.DTTF_phi1[L1Mu_DTTF_index] != 99
has_DT_MB2 = treeHits.DTTF_phi2[L1Mu_DTTF_index] != 99
has_DT_MB3 = treeHits.DTTF_phi3[L1Mu_DTTF_index] != 99
has_DT_MB4 = treeHits.DTTF_phi4[L1Mu_DTTF_index] != 99
if hasMB1Cut and not has_DT_MB1: continue
if hasMB2Cut and not has_DT_MB2: continue
if hasMB3Cut and not has_DT_MB3: continue
if hasMB4Cut and not has_DT_MB4: continue
nDTStubs = has_DT_MB1 + has_DT_MB2 + has_DT_MB3 + has_DT_MB4
if nDTStubs < stubCut: continue
#print "is dt muon"
#print "\t", nDTStubs
## define the L1Mu objects!!
is_CSC_Muon = (1.2 <= abs(L1Mu_eta) and abs(L1Mu_eta) <= 2.4) and L1Mu_CSCTF_index != -1
is_DT_Muon = (0.0 <= abs(L1Mu_eta) and abs(L1Mu_eta) <= 0.9) and L1Mu_DTTF_index != -1
is_DTCSC_Muon = (0.9 <= abs(L1Mu_eta) and abs(L1Mu_eta) <= 1.2) and (L1Mu_DTTF_index != -1 or L1Mu_CSCTF_index != -1)
## muon must be DT or CSC (no RPC)
#if not (is_DT_Muon or is_CSC_Muon):
# continue
if False:
print "\t\tMatched: L1Mu", "pt", L1Mu_pt, "eta", L1Mu_eta,
print "phi", L1Mu_phi, "Quality", L1Mu_quality, "L1mu_bx", L1Mu_bx,
print "L1Mu_CSCTF_index", L1Mu_CSCTF_index, "nCSCStubs", nCSCStubs,
print "L1Mu_DTTF_index", L1Mu_DTTF_index, "nDTStubs", nDTStubs,
#print "Displaced_L1Mu_pt", Displaced_L1Mu_pt
print
#if L1Mu_bx==0:
# print k,i, "pt", L1Mu_pt, "eta", L1Mu_eta, "phi", L1Mu_phi, "bx", L1Mu_bx, "quality", L1Mu_quality
## calculate the max pT for the muons that pass the criteria
if L1Mu_pt > max_prompt_L1Mu_pt:
#print "old L1Mu pt ",max_prompt_L1Mu_pt, " new ",L1Mu_pt
iL1Mu = i
max_prompt_L1Mu_pt = L1Mu_pt
max_prompt_L1Mu_eta = L1Mu_eta
return max_prompt_L1Mu_pt, max_prompt_L1Mu_eta, iL1Mu
def fillDPhiHistogram_MaxPromptPt( h_dphi_ME11_ME21, treeHits, etaCutMin, etaCutMax, npar, timesCharge = True, min_pt = 0, max_pt = 999,
stubCut=2,
qualityCut=4,
hasME1Cut=False,
hasME2Cut=False):
## check if this event has L1Mus
prompt_L1Mu_pt, prompt_L1Mu_eta, i = getMaxPromptPtEtaL1MuEvent(treeHits,
True,
etaCutMin,
etaCutMax,
stubCut,
qualityCut,
hasME1Cut,
hasME2Cut)
if len(list(treeHits.L1Mu_pt))<i or i<0:
#print "len ",len(list(treeHits.L1Mu_pt))," selectedL1Mu ",i
return
#else:
#print "len ",len(list(treeHits.L1Mu_pt))," selectedL1Mu ",i," Prompt_pt ",prompt_L1Mu_pt," eta ", prompt_L1Mu_eta," histo ",h_dphi_ME11_ME21.GetName()
L1Mu_pt = treeHits.L1Mu_pt[i]
L1Mu_eta = treeHits.L1Mu_eta[i]
L1Mu_phi = treeHits.L1Mu_phi[i]
L1Mu_bx = treeHits.L1Mu_bx[i]
L1Mu_quality = treeHits.L1Mu_quality[i]
L1Mu_CSCTF_index = treeHits.L1Mu_CSCTF_index[i]
L1Mu_DTTF_index = treeHits.L1Mu_DTTF_index[i]
L1Mu_charge = treeHits.L1Mu_charge[i]
#print "L1Mu_charge", L1Mu_charge
if L1Mu_pt < min_pt or L1Mu_pt > max_pt: return
## eta cut
if not (etaCutMin <= abs(L1Mu_eta) and abs(L1Mu_eta) <= etaCutMax): return
## quality cut
if L1Mu_quality < qualityCut: return
## BX cut
if abs(L1Mu_bx) != 0 and True: return
## not a CSC muon
if L1Mu_CSCTF_index == -1: return
## parity cut
CSCTF_ch1 = treeHits.CSCTF_ch1[L1Mu_CSCTF_index]
CSCTF_ch2 = treeHits.CSCTF_ch2[L1Mu_CSCTF_index]
CSCTF_isEven1 = CSCTF_ch1%2==0
CSCTF_isEven2 = CSCTF_ch2%2==0
#print "npar required ",npar, " get_parity ",get_parity_ME11_ME21(CSCTF_isEven1, CSCTF_isEven2)
if get_parity_ME11_ME21(CSCTF_isEven1, CSCTF_isEven2) != npar: return
GE11_phi_L1 = treeHits.GE11_phi_L1[L1Mu_CSCTF_index]
GE11_phi_L2 = treeHits.GE11_phi_L2[L1Mu_CSCTF_index]
GE21_phi_L1 = treeHits.GE21_pad2_phi_L1[L1Mu_CSCTF_index]
GE21_phi_L2 = treeHits.GE21_pad2_phi_L2[L1Mu_CSCTF_index]
CSCTF_fit_phi1 = treeHits.CSCTF_fit_phi1[L1Mu_CSCTF_index]
CSCTF_fit_phi2 = treeHits.CSCTF_fit_phi2[L1Mu_CSCTF_index]
GE11_phi = getBestValue(GE11_phi_L1, GE11_phi_L2)
GE21_phi = getBestValue(GE21_phi_L1, GE21_phi_L2)
#GE11_ME11_dphi = - GE11_phi + CSCTF_fit_phi1
#GE21_ME21_dphi = - GE21_phi + CSCTF_fit_phi2
if CSCTF_fit_phi1 == 99 or CSCTF_fit_phi2 == 99:
return
if GE11_phi == 99.:
GE11_ME11_dphi = 99#0.025
elif (timesCharge):
GE11_ME11_dphi = L1Mu_charge*deltaPhi(CSCTF_fit_phi1, GE11_phi)
else:
GE11_ME11_dphi = deltaPhi(CSCTF_fit_phi1, GE11_phi)
if GE21_phi == 99.:
GE21_ME21_dphi = 99#0.025
elif (timesCharge):
GE21_ME21_dphi = L1Mu_charge*deltaPhi(CSCTF_fit_phi2, GE21_phi)
else :
GE21_ME21_dphi = deltaPhi(CSCTF_fit_phi2, GE21_phi)
GE11_dphi_v2 = treeHits.CSCTF_gemdphi1[L1Mu_CSCTF_index]
GE21_dphi_v2 = treeHits.CSCTF_gemdphi2[L1Mu_CSCTF_index]
DPhi = abs(treeHits.CSCTF_L1_DPhi12_GE21[L1Mu_CSCTF_index])
#if ((npar==3 or npar ==1) and npar == get_parity_ME11_ME21(CSCTF_isEven1, CSCTF_isEven2)):
#print "MaxPromptPt, npar required ",npar, " get_parity ",get_parity_ME11_ME21(CSCTF_isEven1, CSCTF_isEven2)," GE11_phi ", GE11_phi," CSCTF_fit_phi1 ",CSCTF_fit_phi1," GE21_phi ", GE21_phi, " CSCTF_fit_phi2 ",CSCTF_fit_phi2," dphi_st2 ",GE21_ME21_dphi
#print "another verioin treeHits.CSCTF_gemdphi1 ",GE11_dphi_v2," gemdphi2 ",GE21_dphi_v2," DPhi_dir ",DPhi,"\n"
#print "fill event, dphi_St1 ", GE11_ME11_dphi," dphi_st2 ",GE21_ME21_dphi
h_dphi_ME11_ME21.Fill(GE11_ME11_dphi, GE21_ME21_dphi)
"""
## calculate the displaced L1Mu pT
doComparatorFit = True
Displaced_L1Mu_pt = pt_endcap_position_based_algorithm(treeHits, i, L1Mu_CSCTF_index, doComparatorFit)
isIsolated = is_L1Mu_isolated(treeHits, i, 0.4, 4, 0.12, 0)
"""
##_________________________________________________
def fillDDYAndDPhiHistogram(h_dYY_dphi, treeHits, etaCutMin, etaCutMax, npar, withGE21, timesCharge = True, min_pt = 0, max_pt = 999,
stubCut =2 ,
qualityCut =4,
hasME1Cut = True,
hasME2Cut = True,
hasME3Cut = True,
hasME4Cut = False):
## check if this event has L1Mus
if len(list(treeHits.L1Mu_pt))==0:
return
pts = list(treeHits.L1Mu_pt)
#print "\tN L1Mus event", len(pts)
for i in range(0,len(pts)):
L1Mu_pt = treeHits.L1Mu_pt[i]
L1Mu_eta = treeHits.L1Mu_eta[i]
L1Mu_phi = treeHits.L1Mu_phi[i]
L1Mu_bx = treeHits.L1Mu_bx[i]
L1Mu_quality = treeHits.L1Mu_quality[i]
L1Mu_CSCTF_index = treeHits.L1Mu_CSCTF_index[i]
L1Mu_DTTF_index = treeHits.L1Mu_DTTF_index[i]
L1Mu_charge = treeHits.L1Mu_charge[i]
#print "L1Mu_charge", L1Mu_charge
if L1Mu_pt < min_pt or L1Mu_pt > max_pt: continue
## eta cut
if not (etaCutMin <= abs(L1Mu_eta) and abs(L1Mu_eta) <= etaCutMax): continue
## quality cut
if L1Mu_quality < qualityCut: continue
## BX cut
if abs(L1Mu_bx) != 0 and True: continue
## not a CSC muon
if L1Mu_CSCTF_index == -1: continue
has_CSC_ME1 = treeHits.CSCTF_phi1[L1Mu_CSCTF_index] != 99
has_CSC_ME2 = treeHits.CSCTF_phi2[L1Mu_CSCTF_index] != 99
has_CSC_ME3 = treeHits.CSCTF_phi3[L1Mu_CSCTF_index] != 99
has_CSC_ME4 = treeHits.CSCTF_phi4[L1Mu_CSCTF_index] != 99
#if hasME1Cut and not has_CSC_ME1: continue
#if hasME2Cut and not has_CSC_ME2: continue
#if hasME3Cut and not has_CSC_ME3: continue
#if hasME4Cut and not has_CSC_ME4: continue
#ok_GE11 = abs(treeHits.CSCTF_gemdphi1[L1Mu_CSCTF_index])< 99
#ok_GE21 = abs(treeHits.CSCTF_gemdphi2[L1Mu_CSCTF_index])< 99
#if ((1.6 <= abs(L1Mu_eta) and abs(L1Mu_eta) <= 2.15 and not(ok_GE11)) or not(ok_GE21)):
# print "GEMHits failed , CSCTF_gemdphi1 ",treeHits.CSCTF_gemdphi1[L1Mu_CSCTF_index]," CSCTF_gemdphi2: ",treeHits.CSCTF_gemdphi2[L1Mu_CSCTF_index]
#else:
# print "GOOD GEMHits"
#if (1.6 <= abs(L1Mu_eta) and abs(L1Mu_eta) <= 2.15 and not(ok_GE11)): continue
#if (withGE21 and not(ok_GE21)): continue
## parity cut
CSCTF_ch1 = treeHits.CSCTF_ch1[L1Mu_CSCTF_index]
CSCTF_ch2 = treeHits.CSCTF_ch2[L1Mu_CSCTF_index]
CSCTF_ch3 = treeHits.CSCTF_ch3[L1Mu_CSCTF_index]
CSCTF_isEven1 = CSCTF_ch1%2==0
CSCTF_isEven2 = CSCTF_ch2%2==0
CSCTF_isEven3 = CSCTF_ch3%2==0
if get_parity(CSCTF_isEven1, CSCTF_isEven2, CSCTF_isEven3) != npar: continue
#with GE21 or not
DPhi_L1 = 99
DDY123_L1 = 99
if (withGE21):
DPhi_L1 = treeHits.CSCTF_L1_DPhi12_GE21[L1Mu_CSCTF_index]
else:
DPhi_L1 = treeHits.CSCTF_L1_DPhi12_noGE21[L1Mu_CSCTF_index]
DDY123_L1 = treeHits.CSCTF_L1_DDY123[L1Mu_CSCTF_index]
if DDY123_L1 == 99:
continue
if (not (has_CSC_ME1 and has_CSC_ME2 and has_CSC_ME3)):
print "AllTracks, DDY123_L1 ",DDY123_L1," CSCTF_phi1 ",treeHits.CSCTF_phi1[L1Mu_CSCTF_index]," CSCTF_phi2 ",treeHits.CSCTF_phi2[L1Mu_CSCTF_index]," CSCTF_phi3 ",treeHits.CSCTF_phi3[L1Mu_CSCTF_index], " CSCTF_eta1 ",treeHits.CSCTF_eta1[L1Mu_CSCTF_index]," CSCTF_eta2 ",treeHits.CSCTF_eta2[L1Mu_CSCTF_index]," CSCTF_eta3 ",treeHits.CSCTF_eta3[L1Mu_CSCTF_index]
if (timesCharge):
DPhi_L1 = L1Mu_charge*DPhi_L1
DDY123_L1 = L1Mu_charge*DDY123_L1
h_dYY_dphi.Fill(DDY123_L1, DPhi_L1)
def fillDisplacedPtEtaHistogram(ptHistogram,
etaHistogram,
treeHits,
doBXCut,
etaCutMin,
etaCutMax,
stubCut,
qualityCut,
hasMB1Cut=False,
hasMB2Cut=False,
hasMB3Cut=False,
hasMB4Cut=False,
doPositionBased=False,
doDirectionBasedNoGE21=False,
doDirectionBasedGE21=False,
doHybridBasedNoGE21=False,
doHybridBasedGE21=False,
doVeto=False,
vetoType=1):
displaced_L1Mu_pt, displaced_L1Mu_eta, iL1Mu= getMaxDisplacedPtEtaL1MuEvent(treeHits,
doBXCut,
etaCutMin,
etaCutMax,
stubCut,
qualityCut,
hasMB1Cut,
hasMB2Cut,
hasMB3Cut,
hasMB4Cut,
doPositionBased,
doDirectionBasedNoGE21,
doDirectionBasedGE21,
doHybridBasedNoGE21,
doHybridBasedGE21,
doVeto,
vetoType)
#print "check pT ok", displaced_L1Mu_pt
if (displaced_L1Mu_pt>0):
ptHistogram.Fill(displaced_L1Mu_pt)
#print "Max displaced pT event", displaced_L1Mu_pt
if (displaced_L1Mu_pt>=10):
etaHistogram.Fill(abs(displaced_L1Mu_eta))
#print "Max displaced pT eta event", displaced_L1Mu_pt, abs(displaced_L1Mu_eta)
## 1: postion based endcap
## 2: direction based endcap
## 3: hybrid endcap
## 4: direction based barrel
displaced_pt_methods = [1, 2, 3]
def getMaxDisplacedPtEtaL1MuEvent(treeHits,
doBXCut,
etaCutMin,
etaCutMax,
stubCut,
qualityCut,
hasMB1Cut=False,
hasMB2Cut=False,
hasMB3Cut=False,
hasMB4Cut=False,
doPositionBased=False,
doDirectionBasedNoGE21=False,
doDirectionBasedGE21=False,
doHybridBasedNoGE21=False,
doHybridBasedGE21=False,
doVeto=False,
vetoType =1):
max_displaced_L1Mu_pt = 0.0
max_displaced_L1Mu_eta = -99
iL1Mu = -1
minDDY123 = 99
minDPhi12_GE21 = 99
minDPhi12_noGE21 = 99
## check if this event has L1Mus
if len(list(treeHits.L1Mu_pt))==0:
return max_displaced_L1Mu_pt, max_displaced_L1Mu_eta, iL1Mu
pts = list(treeHits.L1Mu_pt)
for i in range(0,len(pts)):
L1Mu_eta = treeHits.L1Mu_eta[i]
L1Mu_bx = treeHits.L1Mu_bx[i]
L1Mu_quality = treeHits.L1Mu_quality[i]
## quality cut
if L1Mu_quality < qualityCut: continue
## BX cut
if abs(L1Mu_bx) != 0 and doBXCut: continue
## check if muon is isolated
if doVeto and (not is_L1Mu_isolated(treeHits, i, vetoType)): continue
## eta cut
#if not (etaCutMin <= abs(L1Mu_eta) and abs(L1Mu_eta) <= etaCutMax): continue
L1Mu_CSCTF_index = treeHits.L1Mu_CSCTF_index[i]
L1Mu_DTTF_index = treeHits.L1Mu_DTTF_index[i]
L1Mu_eta_ME2 = -99
## define the L1Mu objects!!
is_CSC_Muon = (1.2 <= abs(L1Mu_eta) and abs(L1Mu_eta) <= 2.4) and L1Mu_CSCTF_index != -1
is_DT_Muon = (0.0 <= abs(L1Mu_eta) and abs(L1Mu_eta) <= 0.9) and L1Mu_DTTF_index != -1
is_DTCSC_Muon = (0.9 <= abs(L1Mu_eta) and abs(L1Mu_eta) <= 1.2) and (L1Mu_DTTF_index != -1 or L1Mu_CSCTF_index != -1)
## CSC quantities
has_CSC_ME1 = False
has_CSC_ME2 = False
has_CSC_ME3 = False
has_CSC_ME4 = False
has_actual_CSC_ME1 = False
has_actual_CSC_ME2 = False
has_CSC_ME11 = False
has_CSC_ME21 = False
is_CSC_ME11_disabled = False
is_CSC_ME21_disabled = False
GE11_dPhi = 99
CSC_ME1_ch = -1
CSC_ME2_ch = -2
nCSCStubs = 0
CSCTF_eta2 = -99
## DT quantities
has_DT_MB1 = False
has_DT_MB2 = False
has_DT_MB3 = False
has_DT_MB4 = False
nDTStubs = 0
DisplacedL1Mu_pt, DisplacedL1Mu_eta = -1, -1
DDY123 = 99
DPhi12_GE21 = 99
DPhi12_noGE21 = 99
#print L1Mu_CSCTF_index
if L1Mu_CSCTF_index != -1 and L1Mu_CSCTF_index < len(treeHits.CSCTF_phi1):
has_CSC_ME1 = treeHits.CSCTF_phi1[L1Mu_CSCTF_index] != 99
has_CSC_ME2 = treeHits.CSCTF_phi2[L1Mu_CSCTF_index] != 99
has_CSC_ME3 = treeHits.CSCTF_phi3[L1Mu_CSCTF_index] != 99
has_CSC_ME4 = treeHits.CSCTF_phi4[L1Mu_CSCTF_index] != 99
L1Mu_eta_ME2 = treeHits.CSCTF_eta2[L1Mu_CSCTF_index]
has_CSC_ME11 = treeHits.CSCTF_ri1[L1Mu_CSCTF_index] == 1 or treeHits.CSCTF_ri1[L1Mu_CSCTF_index] == 4
has_CSC_ME21 = treeHits.CSCTF_ri2[L1Mu_CSCTF_index] == 1
CSC_ME1_ch = treeHits.CSCTF_ch1[L1Mu_CSCTF_index]
CSC_ME2_ch = treeHits.CSCTF_ch2[L1Mu_CSCTF_index]
CSCTF_eta2 = treeHits.CSCTF_eta2[L1Mu_CSCTF_index]
## check if stubs in station 1 and 2 can really be counted! -- they may be disabled
if not has_CSC_ME11: has_actual_CSC_ME1 = has_CSC_ME1
else: has_actual_CSC_ME1 = (not is_CSC_ME11_disabled)
if not has_CSC_ME21: has_actual_CSC_ME2 = has_CSC_ME2
else: has_actual_CSC_ME2 = (not is_CSC_ME21_disabled)
L1Mu_eta = treeHits.CSCTF_L1_eta_st2[L1Mu_CSCTF_index]
DisplacedL1Mu_pt = -1 #not used
DisplacedL1Mu_eta = L1Mu_eta
DDY123 = abs(treeHits.CSCTF_L1_DDY123[L1Mu_CSCTF_index])
DPhi12_GE21 = abs(treeHits.CSCTF_L1_DPhi12_GE21[L1Mu_CSCTF_index])
DPhi12_noGE21 = abs(treeHits.CSCTF_L1_DPhi12_noGE21[L1Mu_CSCTF_index])
if not (etaCutMin <= abs(L1Mu_eta) and abs(L1Mu_eta) <= etaCutMax): continue
#if doPositionBased: DisplacedL1Mu_pt, DisplacedL1Mu_eta = pt_endcap_position_based_algorithm(treeHits, i, True)
#if doDirectionBasedNoGE21: DisplacedL1Mu_pt, DisplacedL1Mu_eta = pt_endcap_direction_based_algorithm(treeHits, i, False)
#if doDirectionBasedGE21: DisplacedL1Mu_pt, DisplacedL1Mu_eta = pt_endcap_direction_based_algorithm(treeHits, i, True)
#if doHybridBasedNoGE21: DisplacedL1Mu_pt, DisplacedL1Mu_eta = pt_endcap_hybrid_algorithm(treeHits, i, False)
#if doHybridBasedGE21: DisplacedL1Mu_pt, DisplacedL1Mu_eta = pt_endcap_hybrid_algorithm(treeHits, i, True)
#print L1Mu_DTTF_index
if is_DT_Muon:
has_DT_MB1 = treeHits.DTTF_phi1[L1Mu_DTTF_index] != 99 and treeHits.DTTF_phib1[L1Mu_DTTF_index] != 99
has_DT_MB2 = treeHits.DTTF_phi2[L1Mu_DTTF_index] != 99 and treeHits.DTTF_phib2[L1Mu_DTTF_index] != 99
has_DT_MB3 = treeHits.DTTF_phi3[L1Mu_DTTF_index] != 99 and treeHits.DTTF_phib3[L1Mu_DTTF_index] != 99
has_DT_MB4 = treeHits.DTTF_phi4[L1Mu_DTTF_index] != 99 and treeHits.DTTF_phib4[L1Mu_DTTF_index] != 99
if hasMB1Cut and not has_DT_MB1: continue
if hasMB2Cut and not has_DT_MB2: continue
if hasMB3Cut and not has_DT_MB3: continue
if hasMB4Cut and not has_DT_MB4: continue
DisplacedL1Mu_pt, DisplacedL1Mu_eta = -1, -1
## calculate the max pT for the muons that pass the criteria
#if DisplacedL1Mu_pt > max_displaced_L1Mu_pt:
if (doPositionBased and abs(DDY123) < minDDY123) or (doDirectionBasedNoGE21 and abs(DPhi12_noGE21) < minDPhi12_noGE21) or (doDirectionBasedGE21 and abs(DPhi12_GE21) < minDPhi12_GE21):
iL1Mu = i
max_displaced_L1Mu_pt = DisplacedL1Mu_pt
max_displaced_L1Mu_eta = DisplacedL1Mu_eta
if abs(DDY123) < minDDY123: minDDY123 = abs(DDY123)
if abs(DPhi12_noGE21) < minDPhi12_noGE21: minDPhi12_noGE21 = abs(DPhi12_noGE21)
if abs(DPhi12_GE21) < minDPhi12_GE21: minDPhi12_GE21 = abs(DPhi12_GE21)
return max_displaced_L1Mu_pt, max_displaced_L1Mu_eta, iL1Mu
##_________________________________________________
def fillDDYAndDPhiHistogram_MaxDisplacedPt( h_dYY_dphi, treeHits, etaCutMin, etaCutMax, npar, withGE21, timesCharge = True, min_pt = 0, max_pt = 999,
stubCut =2 ,
qualityCut =4,
hasME1Cut = True,
hasME2Cut = True,
hasME3Cut = True,
hasME4Cut = False):
displaced_L1Mu_pt, displaced_L1Mu_eta, i = getMaxDisplacedPtEtaL1MuEvent(treeHits,
True,
etaCutMin,
etaCutMax,
3,
4,
False,
False,
False,
False,
True,
False,
False,#dodirectionbasewithGE21
False,
False,
False,
1)
## check if this event has L1Mus
if len(list(treeHits.L1Mu_pt))<i or i<0:
#print "len ",len(list(treeHits.L1Mu_pt))," selectedL1Mu ",i," list of pt ",list(treeHits.L1Mu_pt)
return
#else:
#print "len ",len(list(treeHits.L1Mu_pt))," selectedL1Mu ",i," displaced_pt ", displaced_L1Mu_pt," eta ", displaced_L1Mu_eta," histo ",h_dYY_dphi.GetName()
L1Mu_pt = treeHits.L1Mu_pt[i]
L1Mu_eta = treeHits.L1Mu_eta[i]
L1Mu_phi = treeHits.L1Mu_phi[i]
L1Mu_bx = treeHits.L1Mu_bx[i]
L1Mu_quality = treeHits.L1Mu_quality[i]
L1Mu_CSCTF_index = treeHits.L1Mu_CSCTF_index[i]
L1Mu_DTTF_index = treeHits.L1Mu_DTTF_index[i]
L1Mu_charge = treeHits.L1Mu_charge[i]
#print "L1Mu_charge", L1Mu_charge
if L1Mu_pt< min_pt or L1Mu_pt > max_pt: return
## eta cut
if not (etaCutMin <= abs(L1Mu_eta) and abs(L1Mu_eta) <= etaCutMax): return
## quality cut
if L1Mu_quality < qualityCut: return
## BX cut
if abs(L1Mu_bx) != 0 and True: return
## not a CSC muon
if L1Mu_CSCTF_index == -1: return
has_CSC_ME1 = treeHits.CSCTF_phi1[L1Mu_CSCTF_index] != 99
has_CSC_ME2 = treeHits.CSCTF_phi2[L1Mu_CSCTF_index] != 99
has_CSC_ME3 = treeHits.CSCTF_phi3[L1Mu_CSCTF_index] != 99
has_CSC_ME4 = treeHits.CSCTF_phi4[L1Mu_CSCTF_index] != 99
#if hasME1Cut and not has_CSC_ME1: return
#if hasME2Cut and not has_CSC_ME2: return
#if hasME3Cut and not has_CSC_ME3: return
#if hasME4Cut and not has_CSC_ME4: return
#ok_GE11 = abs(treeHits.CSCTF_gemdphi1[L1Mu_CSCTF_index])< 99
#ok_GE21 = abs(treeHits.CSCTF_gemdphi2[L1Mu_CSCTF_index])< 99
#if (1.6 <= abs(L1Mu_eta) and abs(L1Mu_eta) <= 2.15 and not(ok_GE11)): return
#if (withGE21 and not(ok_GE21)): return
## parity cut
CSCTF_ch1 = treeHits.CSCTF_ch1[L1Mu_CSCTF_index]
CSCTF_ch2 = treeHits.CSCTF_ch2[L1Mu_CSCTF_index]
CSCTF_ch3 = treeHits.CSCTF_ch3[L1Mu_CSCTF_index]
CSCTF_isEven1 = CSCTF_ch1%2==0
CSCTF_isEven2 = CSCTF_ch2%2==0
CSCTF_isEven3 = CSCTF_ch3%2==0
if get_parity(CSCTF_isEven1, CSCTF_isEven2, CSCTF_isEven3) != npar: return
#with GE21 or not
DPhi_L1 = 99
DDY123_L1 = 99
if (withGE21):
DPhi_L1 = treeHits.CSCTF_L1_DPhi12_GE21[L1Mu_CSCTF_index]
else:
DPhi_L1 = treeHits.CSCTF_L1_DPhi12_noGE21[L1Mu_CSCTF_index]
DDY123_L1 = treeHits.CSCTF_L1_DDY123[L1Mu_CSCTF_index]
if DDY123_L1 == 99: return
if (not (has_CSC_ME1 and has_CSC_ME2 and has_CSC_ME3)):
print "MaxDisplacedPt DDY123_L1 ",DDY123_L1," CSCTF_phi1 ",treeHits.CSCTF_phi1[L1Mu_CSCTF_index]," CSCTF_phi2 ",treeHits.CSCTF_phi2[L1Mu_CSCTF_index]," CSCTF_phi3 ",treeHits.CSCTF_phi3[L1Mu_CSCTF_index], " CSCTF_eta1 ",treeHits.CSCTF_eta1[L1Mu_CSCTF_index]," CSCTF_eta2 ",treeHits.CSCTF_eta2[L1Mu_CSCTF_index]," CSCTF_eta3 ",treeHits.CSCTF_eta3[L1Mu_CSCTF_index]
if (timesCharge):
DPhi_L1 = L1Mu_charge*DPhi_L1
DDY123_L1 = L1Mu_charge*DDY123_L1
#print "fillevent DDY123_L1 ",DDY123_L1," DPhi_L1 ",DPhi_L1
h_dYY_dphi.Fill(DDY123_L1, DPhi_L1)
## match the L1Mu to EMTF Mu
def get_L1Mu_CSCTF_index(SIM_L1Mu_index, treeHits):
if SIM_L1Mu_index == 999: return 99
L1Mu_CSCTF_index = -1
maxdR = 999
for jjj in range(0,len(treeHits.CSCTF_pt)):
if treeHits.CSCTF_bx[jjj] != 0:
continue
#print "\t", treeHits.CSCTF_pt[jjj], treeHits.CSCTF_eta[jjj], treeHits.CSCTF_phi[jjj], treeHits.CSCTF_bx[jjj]
dEta = treeHits.CSCTF_eta[jjj] - treeHits.L1Mu_eta[SIM_L1Mu_index]
dPhi = treeHits.CSCTF_phi[jjj] - treeHits.L1Mu_phi[SIM_L1Mu_index]
dR = sqrt(dEta * dEta + dPhi*dPhi)
if dR < maxdR and dR < 0.1:
maxdR = dR
L1Mu_CSCTF_index = jjj
return L1Mu_CSCTF_index
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment