Skip to content

Instantly share code, notes, and snippets.

@tahuang1991

tahuang1991/makeL1HybridTriggerRatePlots_tree.py

Created February 23, 2017 20:50
Show Gist options
  • Save tahuang1991/baefc4c678c85be35ed0c40ab9ab2169 to your computer and use it in GitHub Desktop.
Save tahuang1991/baefc4c678c85be35ed0c40ab9ab2169 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
makeL1HybridTriggerRatePlots_tree.py
# Run quiet mode
import sys
sys.argv.append( '-b' )
import ROOT
ROOT.gROOT.SetBatch(1)
from Helpers import *
from hybridAlgorithmPtAssignment import *
from TTTrackIsolation import *
ROOT.gErrorIgnoreLevel=1001
from ROOT import *
import random
import numpy as n
#______________________________________________________________________________
if __name__ == "__main__":
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170105_v2"; pu = 'PU140'; eff = False
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170110"; pu = 'PU140'; eff = False
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170119"; pu = 'PU140'; eff = False
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170123"; pu = 'PU140'; eff = False
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170124"; pu = 'PU140'; eff = False
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170125"; pu = 'PU140'; eff = False
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170126"; pu = 'PU140'; eff = False
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170131"; pu = 'PU140'; eff = False
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170207"; pu = 'PU140'; eff = False
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170207_v2"; pu = 'PU140'; eff = False
label = "Neutrino_Pt2to20_gun_TTI2023Upg14D_PU140bx25_ILT_SLHC14_20170207_v2_1k"; pu = 'PU140'; eff = False
label = "MinBias_L1_CMSSW90X_10M"; pu = 'PU0'; eff = False
if len(sys.argv)>=2:
label = sys.argv[1]
## extension for figures - add more?
ext = ".png"
## Style
#gStyle.SetStatStyle(0)
print "Making the plots"
set_style()
doTest = True
ch = TChain("DisplacedL1MuFilter_PhaseIIGE21/L1MuTree")
location0 = '/eos/uscms/store/user/lpcgem/MinBias_TuneCUETP8M1_14TeV-pythia8/MinBias_TuneCUETP8M1_14TeV_ANA_20170210_v1/170211_061226/0000/'
location0 = '/eos/uscms/store/user/tahuang/CMSSW620SLHC/RateSample/'
#treeHits = addfiles(ch, dirname=location0, ext=".root")
##treeHits = addfiles(ch, "out_ana_rate.root")
#treeHits = addfiles(ch, "out_ana_0.root")
#treeHits = addfiles(ch, "out_ana_1.root")
treeHits = addfiles(ch, "out_ana_2.root")
#treeHits = addfiles(ch, location0)
#treeHits = addfiles(ch, "ratesample/")
targetDir = label + "/"
if not os.path.exists(targetDir):
os.makedirs(targetDir)
verbose = False
## copy index file
#import shutil
#shutil.copy2('index.php', targetDir + 'index.php')
def displacedL1MuHybridTriggerRate():
## helper functions to make many plots!!
mapTH1F = ROOT.std.map("string,TH1F")()
mapTH2F = ROOT.std.map("string,TH2F")()
def addPlotToMapTH1F(name,nBin,minBin,maxBin):
mapTH1F[name] = TH1F(name,"",nBin,minBin,maxBin)
def addPlotToMapTH1F_v2(name,bins):
mapTH1F[name] = TH1F(name,"",len(bins)-1, bins)
def addManyPlotsToTH1F(ptbins, etabins, *arg):
for i in range(1,len(arg)):
addPlotToMapTH1F_v2(arg[i].replace("rate_", "rate_pt_"), ptbins)
addPlotToMapTH1F_v2(arg[i].replace("rate_", "rate_eta_"), etabins)
## binning
ptbin = [
1, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0,
10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0,
45.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 120.0, 140.0]
myptbin = np.asarray(ptbin)
nmyptbin = len(myptbin) - 1
etabin = [
0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
1.0, 1.05, 1.1, 1.15, 1.2, 1.25, 1.3, 1.35, 1.4, 1.45, 1.5, 1.55, 1.6, 1.65, 1.7, 1.75, 1.8, 1.85, 1.9, 1.95,
2.0, 2.05, 2.1, 2.15, 2.2, 2.25, 2.3, 2.35, 2.4, 2.45, 2.5]
myetabin = np.asarray(etabin)
L1MuRatetree = TTree("L1MuTriggerRate", "L1MuTriggerRate")
position_pt = n.zeros(1, dtype=float)
direction_pt_noGE21 = n.zeros(1, dtype=float)
direction_pt_GE21 = n.zeros(1, dtype=float)
hybrid_pt_noGE21 = n.zeros(1, dtype=float)
hybrid_pt_GE21 = n.zeros(1, dtype=float)
dphi_dir12_noGE21 = n.zeros(1, dtype=float)
dphi_dir12_GE21 = n.zeros(1, dtype=float)
ddY123 = n.zeros(1, dtype=float)
gemdphi_st1 = n.zeros(1, dtype=float)
gemdphi_st2 = n.zeros(1, dtype=float)
npar= n.zeros(1, dtype=int)
npar_bending= n.zeros(1, dtype=int)
passGE11 = n.zeros(1, dtype=bool)
passGE21 = n.zeros(1, dtype=bool)
maxPositionPt = n.zeros(1, dtype=bool)
maxDirectionGE21Pt = n.zeros(1, dtype=bool)
maxDirectionNoGE21Pt = n.zeros(1, dtype=bool)
maxPromptPt = n.zeros(1, dtype=bool)
#L1Mu track
L1Mu_pt = n.zeros(1, dtype=float)
L1Mu_charge = n.zeros(1, dtype=int)
L1Mu_eta = n.zeros(1, dtype=float)
L1Mu_eta_st2 = n.zeros(1, dtype=float)
L1Mu_phi = n.zeros(1, dtype=float)
L1Mu_bx = n.zeros(1, dtype=int)
L1Mu_nstubs = n.zeros(1, dtype=int)
L1Mu_quality = n.zeros(1, dtype=int)
L1MuRatetree.Branch("position_pt",position_pt,"position_pt/D")
L1MuRatetree.Branch("direction_pt_GE21",direction_pt_GE21,"direction_pt_GE21/D")
L1MuRatetree.Branch("direction_pt_noGE21",direction_pt_noGE21,"direction_pt_noGE21/D")
L1MuRatetree.Branch("hybrid_pt_GE21",hybrid_pt_GE21,"hybrid_pt_GE21/D")
L1MuRatetree.Branch("hybrid_pt_noGE21",hybrid_pt_noGE21,"hybrid_pt_noGE21/D")
L1MuRatetree.Branch("dphi_dir12_GE21",dphi_dir12_GE21,"dphi_dir12_GE21/D")
L1MuRatetree.Branch("dphi_dir12_noGE21",dphi_dir12_noGE21,"dphi_dir12_noGE21/D")
L1MuRatetree.Branch("ddY123",ddY123,"ddY123/D")
L1MuRatetree.Branch("gemdphi_st1",gemdphi_st1,"gemdphi_st1/D")
L1MuRatetree.Branch("gemdphi_st2",gemdphi_st2,"gemdphi_st2/D")
L1MuRatetree.Branch("npar",npar,"npar/I")
L1MuRatetree.Branch("npar_bending",npar_bending,"npar_bending/I")
L1MuRatetree.Branch("passGE11",passGE11,"passGE11/B")
L1MuRatetree.Branch("passGE21",passGE21,"passGE21/B")
L1MuRatetree.Branch("maxPositionPt",maxPositionPt,"maxPositionPt/B")
L1MuRatetree.Branch("maxDirectionGE21Pt",maxDirectionGE21Pt,"maxDirectionGE21Pt/B")
L1MuRatetree.Branch("maxDirectionNoGE21Pt",maxDirectionNoGE21Pt,"maxDirectionNoGE21Pt/B")
L1MuRatetree.Branch("maxPromptPt",maxPromptPt,"maxPromptPt/B")
L1MuRatetree.Branch("L1Mu_pt",L1Mu_pt,"L1Mu_pt/D")
L1MuRatetree.Branch("L1Mu_charge",L1Mu_charge,"L1Mu_charge/I")
L1MuRatetree.Branch("L1Mu_eta",L1Mu_eta,"L1Mu_eta/D")
L1MuRatetree.Branch("L1Mu_eta_st2",L1Mu_eta_st2,"L1Mu_eta_st2/D")
L1MuRatetree.Branch("L1Mu_phi",L1Mu_phi,"L1Mu_phi/D")
L1MuRatetree.Branch("L1Mu_bx",L1Mu_bx,"L1Mu_bx/I")
L1MuRatetree.Branch("L1Mu_nstubs",L1Mu_nstubs,"L1Mu_nstubs/I")
L1MuRatetree.Branch("L1Mu_quality",L1Mu_quality,"L1Mu_quality/I")
#L1MuRatetree.Branch()
h_eventcount = TH1F("h_eventcount","",10,0,10)
etapartition = [1.2,1.4,1.6,1.8,2.0,2.15]
nparity = 4
maxEntries = ch.GetEntries()
if doTest:
maxEntries = 10000
nEvents3stationPassPrompt = 0
nEvents3stationPassDisplaced = 0
nEvents = maxEntries
print "nEvents", nEvents
for k in range(0,nEvents):
if k%1000==0: print "Processing event", k
h_eventcount.Fill(1)
ch.GetEntry(k)
treeHits = ch
def initbranches():
#init branches
position_pt[0] = -1
direction_pt_noGE21[0] = -1
direction_pt_GE21[0] = -1
hybrid_pt_noGE21[0] = -1
hybrid_pt_GE21[0] = -1
dphi_dir12_noGE21[0] = 99
dphi_dir12_GE21[0] = 99
ddY123[0] = 99
gemdphi_st1[0] = 99
gemdphi_st2[0] = 99
npar[0] = -1
npar_bending[0]= -1
passGE11[0] = False
passGE21[0] = False
maxPromptPt[0] = False
maxPositionPt[0] = False
maxDirectionGE21Pt[0] = False
maxDirectionNoGE21Pt[0] = False
L1Mu_pt[0] = -1
L1Mu_eta[0]= -9
L1Mu_eta_st2[0]= -9
L1Mu_phi[0] = -9
L1Mu_charge[0] = -9
L1Mu_bx[0] = 0
L1Mu_quality[0] = -1
initbranches()
doBXCut = True
etaCutMin = 1.6
etaCutMax = 2.4
stubCut = 2
qualityCut = 4
hasME1Cut = True
hasME2Cut = True
hasME3Cut = True
hasME4Cut = False
hasME11Cut = True
hasME21Cut = True
hasGE11Cut = False
hasGE21Cut = False
hasME11ME21Cut = True
hasGE11GE21Cut = True
prompt_L1Mu_pt, prompt_L1Mu_eta, prompt_iL1Mu = getMaxPromptPtEtaL1MuEvent(treeHits,
doBXCut,
etaCutMin,
etaCutMax,
stubCut,
qualityCut,
hasME1Cut,
hasME2Cut,
hasME3Cut,
hasME4Cut,
hasME11Cut,
hasME21Cut,
hasGE11Cut,
hasGE21Cut,
hasME11ME21Cut,
hasGE11GE21Cut,
0,#ME11FailRate,
0,#ME21FailRate,
False,#hasMB1Cut,
False,#hasMB2Cut,
False,#hasMB3Cut,
False#hasMB4Cut
)
#get max pt from position based
displaced_L1Mu_pt, displaced_L1Mu_eta, displaced_iL1Mu = getMaxDisplacedPtEtaL1MuEvent(treeHits,
doBXCut,
etaCutMin,
etaCutMax,
stubCut,
qualityCut,
False,#hasMB1Cut,
False,#hasMB2Cut,
False,#hasMB3Cut,
False,#hasMB4Cut,
True,#doPositionBased,
True,#doDirectionBasedNoGE21,
True,#doDirectionBasedGE21,
False,#doHybridBasedNoGE21,
False,#doHybridBasedGE21,
False,#doVeto,
False#vetoType
)
displaced_L1Mu_directionNoGE21_pt, displaced_L1Mu_directionNoGE21_eta, displaced_directionNoGE21_iL1Mu = getMaxDisplacedPtEtaL1MuEvent(treeHits,
doBXCut,
etaCutMin,
etaCutMax,
stubCut,
qualityCut,
False,#hasMB1Cut,
False,#hasMB2Cut,
False,#hasMB3Cut,
False,#hasMB4Cut,
#True,#doPositionBased,
True,#doDirectionBasedNoGE21,
False,#doDirectionBasedGE21,
False,#doHybridBasedNoGE21,
False,#doHybridBasedGE21,
False,#doVeto,
False#vetoType
)
# ptHistogram.Fill(prompt_L1Mu_pt)
## get the max value of the momentum
pts = list(treeHits.L1Mu_pt)
## ignore events without L1Mu
if len(pts)==0: continue
minQuality = 4
for i in range(0,len(pts)):
initbranches()
if i==prompt_iL1Mu:
maxPromptPt[0] = True
if i==displaced_iL1Mu:
maxPositionPt[0] = True
if i==displaced_directionNoGE21_iL1Mu:
maxDirectionNoGE21Pt[0] = True
L1Mu_pt[0] = treeHits.L1Mu_pt[i]
L1Mu_eta[0] = treeHits.L1Mu_eta[i]
L1Mu_phi[0] = treeHits.L1Mu_phi[i]
L1Mu_bx[0] = treeHits.L1Mu_bx[i]
L1Mu_quality[0] = treeHits.L1Mu_quality[i]
L1Mu_CSCTF_index = treeHits.L1Mu_CSCTF_index[i]
L1Mu_DTTF_index = treeHits.L1Mu_DTTF_index[i]
L1Mu_charge[0] = treeHits.L1Mu_charge[i]
#L1Mu_nstubs[0] =
## eta cut
if not (etaCutMin <= abs(L1Mu_eta[0]) and abs(L1Mu_eta[0]) <= etaCutMax): continue
## quality cut
if L1Mu_quality[0] < qualityCut: continue
## BX cut
if abs(L1Mu_bx[0]) != 0 and doBXCut: continue
## not a CSC muon
if L1Mu_CSCTF_index == -1: continue
dphi_dir12_GE21[0] = float(treeHits.CSCTF_L1_DPhi12_GE21[L1Mu_CSCTF_index])
dphi_dir12_noGE21[0] = treeHits.CSCTF_L1_DPhi12_noGE21[L1Mu_CSCTF_index]
ddY123[0] = float(treeHits.CSCTF_L1_DDY123[L1Mu_CSCTF_index])
position_pt[0] = treeHits.CSCTF_L1_position_pt[L1Mu_CSCTF_index]
direction_pt_noGE21[0] = treeHits.CSCTF_L1_direction_pt_noGE21[L1Mu_CSCTF_index]
direction_pt_GE21[0] = treeHits.CSCTF_L1_direction_pt_GE21[L1Mu_CSCTF_index]
hybrid_pt_noGE21[0] = treeHits.CSCTF_L1_hybrid_pt_noGE21[L1Mu_CSCTF_index]
hybrid_pt_GE21[0] = treeHits.CSCTF_L1_hybrid_pt_GE21[L1Mu_CSCTF_index]
L1Mu_eta_st2[0] = treeHits.CSCTF_eta2[L1Mu_CSCTF_index]
#double check
#if maxPromptPt[0]:
# print "L1Mu_pt ",L1Mu_pt[0]," L1Mu_eta ", L1Mu_eta[0]," prompt_L1Mu_pt ",prompt_L1Mu_pt," prompt_L1Mu_eta ",prompt_L1Mu_eta
#if maxPositionPt[0]:
# print "displaced_pt ",position_pt[0]," L1Mu_eta ", L1Mu_eta_st2[0]," displaced_L1Mu_pt ",displaced_L1Mu_pt," displaced_L1Mu_eta ",displaced_L1Mu_eta
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_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
npar_bending[0] = get_parity_ME11_ME21(CSCTF_isEven1, CSCTF_isEven2)
npar[0] = get_parity(CSCTF_isEven1, CSCTF_isEven2, CSCTF_isEven3)
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)
gemdphi_st1[0] = deltaPhi(CSCTF_fit_phi1, GE11_phi)
gemdphi_st2[0] = deltaPhi(CSCTF_fit_phi2, GE21_phi)
if CSCTF_fit_phi1 == 99 or GE11_phi == 99:
gemdphi_st1[0] = 99
if CSCTF_fit_phi2 == 99 or GE21_phi == 99:
gemdphi_st2[0] = 99
if gemdphi_st1[0] != 99:
passGE11[0] = passDPhicutTFTrack(1, CSCTF_ch1, gemdphi_st1[0], L1Mu_pt[0])
if gemdphi_st2[0] != 99:
passGE21[0] = passDPhicutTFTrack(2, CSCTF_ch2, gemdphi_st2[0], L1Mu_pt[0])
#which gem dphi to use??
GE11_dphi_v2 = treeHits.CSCTF_gemdphi1[L1Mu_CSCTF_index]
GE21_dphi_v2 = treeHits.CSCTF_gemdphi2[L1Mu_CSCTF_index]
print "gemdphi_st1[0] ",gemdphi_st1[0]," gemdphi_st2 ",gemdphi_st2[0], " ddY123 ",ddY123[0], " dphi_dir_GE21 ",dphi_dir12_GE21[0], " dphi_dir_noGE21 ",dphi_dir12_noGE21[0]
L1MuRatetree.Fill()
## make plots
def makeDPhiPlot(h, xaxis, yaxis):
c = TCanvas("c","c",800,600)
c.Clear()
title = h.GetName()
gStyle.SetOptStat(111111)
gStyle.SetTitleBorderSize(0);
#gStyle.SetPadLeftMargin(0.126);
#gStyle.SetPadRightMargin(0.10);
#gStyle.SetPadTopMargin(0.06);
#gStyle.SetPadBottomMargin(0.13);
gPad.SetTickx(1)
gPad.SetTicky(1)
gPad.SetLogx(0);
gPad.SetLogy(0)
gPad.SetGridx(1);
gPad.SetGridy(1);
gStyle.SetTitleStyle(0)
gStyle.SetTitleAlign(13) ##// coord in top left
gStyle.SetTitleX(0.)
gStyle.SetTitleY(1.)
gStyle.SetTitleW(1)
gStyle.SetTitleH(0.058)
gStyle.SetTitleBorderSize(0)
gStyle.SetPadLeftMargin(0.126)
gStyle.SetPadRightMargin(0.10)
gStyle.SetPadTopMargin(0.06)
gStyle.SetPadBottomMargin(0.13)
gStyle.SetOptStat(0)
gStyle.SetMarkerStyle(1)
h.GetXaxis().SetTitleSize(0.05)
h.GetYaxis().SetTitleSize(0.05)
h.GetXaxis().SetLabelSize(0.05)
h.GetYaxis().SetLabelSize(0.05)
gStyle.SetTitleFontSize(0.065)
h.SetTitle(" #scale[1.4]{#font[61]{CMS}} #font[52]{Simulation preliminary} 14 TeV, 140 PU")
#h.SetStats(0)
#h.GetYaxis().SetTitle("charge*#Delta#Phi(GE1/1-ME1/1)")
#h.GetXaxis().SetTitle("charge*#Delta#Phi(GE2/1-ME2/1)")
h.GetYaxis().SetTitle("%s"%(yaxis))
h.GetXaxis().SetTitle("%s"%(xaxis))
h.Draw("COLZ")
c.SaveAs(targetDir + title + ".png")
#c.SaveAs(targetDir + title + ".pdf")
c.SaveAs(targetDir + title + ".C")
#makeDPhiPlot(h_dphi_ME11_ME21, "dphi_ME11_ME21")
c = TCanvas("c","c",800,600)
c.Clear()
h_eventcount.Draw()
c.SaveAs("h_eventcount.png")
#make plots before write into root
targetroot = TFile("target.root","RECREATE")
h_eventcount.Write()
#h_dphi_ME11_ME21.Write()
L1MuRatetree.Write()
targetroot.Close()
## trigger rate plots
displacedL1MuHybridTriggerRate()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment