L1TCalorimeter/python/studyJets.py

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