L1TCalorimeter/python/studyJets.py

0001 #!/usr/bin/env python
0002
0003 from __future__ import print_function
0004 import os
0005 import sys
0006 import math
0007
0008
0009 ####################################################################################
0010 #  Python script for studying eta dependence of jet sizes and seed centrality      #
0011 #  Author: Aaron Bundock (aaron.*nospamthankyamaam*bundock@cern.ch)                #
0012 #  Date: 15/01/18                                                                  #
0013 #                                                                                  #
0014 ####################################################################################
0015
0016
0017
0018 towerEtaWidths = [0, 0.087, 0.087, 0.087, 0.087, 0.087, 0.087, 0.087, 0.087, 0.087, 0.087,                 # 0-10
0019                      0.087, 0.087, 0.087, 0.087, 0.087, 0.087, 0.087, 0.087, 0.087, 0.087,                 # 11-20
0020                      0.090, 0.100, 0.113, 0.129, 0.150, 0.178, 0.150, 0.350,                               # 21-28
0021                      0.132, 0.175, 0.176, 0.174, 0.176, 0.174, 0.177, 0.173, 0.175, 0.177, 0.173, 0.302]   # 30-41
0022 #0.150, 0.180, 0.170, 0.180, 0.170, 0.180, 0.170, 0.180, 0.180, 0.180, 0.180, 0.290]
0023
0024
0025
0026 def getJetProperties(jetSeed,etaFwd,etaCen):
0027
0028     jetSize = 0
0029     etaFwdSize = 0
0030     etaCenSize = 0
0031
0032     for ring in xrange(jetSeed-etaCen,jetSeed+etaFwd+1):
0033
0034         if ring < 1:
0035             ring = abs(ring-1)
0036
0037         if ring >= len(towerEtaWidths):
0038             break
0039
0040         if ring < jetSeed:
0041             etaFwdSize += towerEtaWidths[ring]
0042         if ring > jetSeed:
0043             etaCenSize += towerEtaWidths[ring]
0044
0045     jetSize = (etaFwdSize + etaCenSize + towerEtaWidths[jetSeed])/0.8
0046     seedCent = etaCenSize/etaFwdSize
0047
0048     jetProps = [jetSize, seedCent]
0049
0050     return jetProps
0051
0052
0053 def printJetProperties(etaRange):
0054
0055     print("Size  \  eta\t", end=' ')
0056     for seedEta in etaRange:
0057         if(seedEta<29):
0058             print(str(seedEta)+"\t\t", end=' ')
0059         else:
0060             print(str(seedEta+1)+"\t\t", end=' ')
0061     print()
0062
0063     for size in jetEtaSizes:
0064
0065         if size%2:
0066             print("   9x"+str(size)+("\t"), end=' ')
0067             for seedEta in etaRange:
0068                 print(("\t"), end=' ')
0069                 etaFwdCen = (size-1)/2
0070                 jetProps = getJetProperties(seedEta,etaFwdCen,etaFwdCen)
0071                 print(("%.2f / %.2f" %(jetProps[0],jetProps[1])), end=' ')
0072             print()
0073
0074         else:
0075             print("   9x"+str(size)+" cen", end=' ')
0076             for seedEta in etaRange:
0077                 print(("\t"), end=' ')
0078                 etaFwd = size/2-1
0079                 etaCen = size/2
0080                 jetProps = getJetProperties(seedEta, etaFwd, etaCen)
0081                 print(("%.2f / %.2f" %(jetProps[0],jetProps[1])), end=' ')
0082             print()
0083             print("   9x"+str(size)+" fwd", end=' ')
0084             for seedEta in etaRange:
0085                 print(("\t"), end=' ')
0086                 etaFwd = size/2
0087                 etaCen = size/2-1
0088                 jetProps = getJetProperties(seedEta, etaFwd, etaCen)
0089                 print(("%.2f / %.2f" %(jetProps[0],jetProps[1])), end=' ')
0090             print()
0091
0092
0093
0094 print("\n")
0095 print("==============================================================================================================================================================================================================================")
0096 print("Eta-dependence of jet sizes & seed centrality")
0097 print("First number  = jet diameter in eta, normalised to 0.8 (offline)")
0098 print("Second number = jet size on forward side of jet seed / jet size on central side of jet seed")
0099 print("Ideally, best jet choice is where both numbers are closest to unity!")
0100 print("9xN cen = larger area on side of jet further from beam pipe")
0101 print("9xN fwd  = larger area on side of jet closer to beam pipe")
0102 print("==============================================================================================================================================================================================================================")
0103
0104
0105 beEtaRange = xrange(16,29)
0106 hfEtaRange = xrange(29,41)
0107
0108
0109 jetEtaSizes = [9,8,7,6,5,4,3]
0110
0111
0112 print("\n BE \n")
0113 printJetProperties(beEtaRange)
0114
0115 print("\n\n HF \n")
0116 printJetProperties(hfEtaRange)
0117
0118 print("\n")
0119 print("==============================================================================================================================================================================================================================")
0120 print("\n\n")