qcd_pythia8_eta_jets.py (raw text file)# Part of =HepSim=: http://atlaswww.hep.anl.gov/hepsim/
# Jython validation script using Jas4pp
# Task: show inclusive pT cross sections for different Eta ranges
# S.Chekanov (ANL)
import sys
from java.awt import Color,Font
from jhplot import HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from java.lang import *
from proto import FileMC # import FileMC
from jhplot.utils import FileList
from hephysics.particle import LParticle
from hepsim import HepSim
from java.util import ArrayList
from hephysics.hepsim import PromcUtil
from hephysics.jet import JetN2
from math import *
ktjet=JetN2(0.5,"antikt",40) # R=0.5, E-mode, anti-KT,min pT=40
ktjet.setDebug(False)
NMax=3 # max number of files per pT. Set to -1 to run over all files
# default location
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/qcd_pythia8_ptall/"
# get data either from local file system or stream them from URL
url=""
TotalEvents=100000 # number of events per sample
if len(sys.argv)>1:
if sys.argv[1].startswith("http"):
flist=HepSim.getList(sys.argv[1])
url=sys.argv[1]
else: flist=FileList.get(sys.argv[1],"promc")
if (len(sys.argv)>2): TotalEvents=int(sys.argv[2])
else:
url=default_www; flist=HepSim.getList(url)
if len(flist)==0: print "Error: No input file!"; sys.exit(0)
else: print "Reading "+str(len(flist))+" files. Nr events= ",TotalEvents
# get list of files with different pT range
xmin=HepSim.getRanges(flist,"_pt","_");
# reduce range
# del xmin[-1]; del xmin[-1];
# define histogram bins in ranges of pT, with 10% shift to avoid biases
bins=[]
for j in xmin:
bins.append(j+j*0.1)
print "Histogram bins",bins
# create file list using pT ranges
filelist=[]
histo1=[] # 0<|eta|<0.5
histo2=[] # 0.5<|eta|<1.0
histo3=[] # 1.0<|eta|<1.5
histo4=[] # 1.5<|eta|<2.0
histo5=[] # 2.0<|eta|<2.5
for i in xmin:
pt=[]
for j in flist:
if (j.find("pt"+str(i)+"_")>-1): pt.append(j)
filelist.append(pt)
histo1.append( H1D("eta05 pt"+str(i),bins) ) # 50 GeV bins
histo2.append( H1D("eta10 pt"+str(i),bins) ) # 50 GeV bins
histo3.append( H1D("eta15 pt"+str(i),bins) ) # 50 GeV bins
histo4.append( H1D("eta20 pt"+str(i),bins) ) # 50 GeV bins
histo5.append( H1D("eta25 pt"+str(i),bins) ) # 50 GeV bins
print "Total number of pT ranges=",len(filelist)
Total=len(filelist)
for ipt in range(Total):
ptlist=filelist[ipt]
print ipt, ") pT range = "+str(xmin[ipt])
cross=0; xcross=0; nev=0; Nfiles=len(ptlist)
if (NMax>-1): Nfiles=NMax # limit to NMax files
for m in range(Nfiles): # loop over all files in this list
file=FileMC(url+ptlist[m]) # get input file
print "Read=",url+ptlist[m]
if (nev>TotalEvents): print "Max Nr of events are done"; break # stop after some events
for i in range(file.size()):
if (nev>TotalEvents): break
nev=nev+1
if (nev%200==0):
if (m==0): print "Event=",nev
else: print "Event=",nev," done=",int((100.0*m)/Nfiles),"%", 'x-cross=',xcross
eve = file.read(i)
pa = eve.getParticles() # particle information
ve = eve.getEvent() # event information
particles=PromcUtil.getParticleDList(file.getHeader(), pa, 1, -1, 1000);
ktjet.buildJets(particles);
# ktjet.printJets();
jets=ktjet.getJetsSorted()
for jet in jets:
pt=jet.perp()
eta=jet.eta()
if (pt>xmin[ipt]):
if (abs(eta)<0.5): histo1[ipt].fill(pt)
if (abs(eta)>=0.5 and abs(eta)<1.0): histo2[ipt].fill(pt)
if (abs(eta)>=1.0 and abs(eta)<1.5): histo3[ipt].fill(pt)
if (abs(eta)>=1.5 and abs(eta)<2.0): histo4[ipt].fill(pt)
if (abs(eta)>=2.0 and abs(eta)<2.5): histo5[ipt].fill(pt)
stat = file.getStat()
xcross=stat.getCrossSectionAccumulated()
cross=cross+xcross;
erro=stat.getCrossSectionErrorAccumulated();
file.close()
cross=cross/Nfiles
lumi=nev/cross # get cross sections
histo1[ipt].scale(1.0/lumi)
histo2[ipt].scale(1.0/lumi)
histo3[ipt].scale(1.0/lumi)
histo4[ipt].scale(1.0/lumi)
histo5[ipt].scale(1.0/lumi)
print ipt, ") Cross section = %.3e pb"%cross, " Lumi=%.3e pb"%lumi
h1=histo1[0] # combine all pT ranges
h1.setTitle("0<|η|<0.5")
h1.setColor(Color.black)
h1.setErrColor(Color.black)
for i in range(1,len(histo1)):
h1=h1.oper(histo1[i],"+")
h2=histo2[0] # combine all pT ranges
h2.setTitle("0.5<|η|<1.0")
for i in range(1,len(histo2)):
h2=h2.oper(histo2[i],"+")
h2.setColor(Color.red)
h2.setErrColor(Color.red)
h3=histo3[0] # combine all pT ranges
h3.setTitle("1.0<|η|<1.5")
for i in range(1,len(histo3)):
h3=h3.oper(histo3[i],"+")
h3.setColor(Color.blue)
h3.setErrColor(Color.blue)
h4=histo4[0] # combine all pT ranges
h4.setTitle("1.5<|η|<2.0")
for i in range(1,len(histo4)):
h4=h4.oper(histo4[i],"+")
h4.setColor(Color.green)
h4.setErrColor(Color.green)
h5=histo5[0] # combine all pT ranges
h5.setTitle("2.0<|η|<2.5")
for i in range(1,len(histo5)):
h5=h5.oper(histo5[i],"+")
h5.setColor(Color.magenta)
h5.setErrColor(Color.magenta)
crossD=h1.sumAllBinHeights() # get total cross section
h1=h1.getDividedByBinWidth() # get differential cross section
c1 = HPlot("HepSim",600,400) # plot histogram
c1.setNameX("p_{T}(jet) [GeV]")
c1.setNameY("d σ / d p_{T} (jet) [pb / GeV]")
c1.visible()
c1.setMarginLeft(100)
c1.setAutoRange()
c1.setRangeX(0,30000)
c1.setRangeY(10e-13,10e+8)
c1.setLogScale(1,True)
c1.draw([h1,h2,h3,h4,h5])
# h1.toTable()
l1=HLabel("Pythia8 QCD σ=%.3e pb"%crossD, 0.5, 0.65, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l1)
l2=HLabel("=HepSim=",0.7,0.84, "NDC")
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
l2.setColor(Color.gray)
c1.add(l2)
# create file/image using name of the file
name="output"
if len(sys.argv[0])>0: name=sys.argv[0].replace(".py","")
file=HBook(name+".jdat","w"); print name+".jdat created"
file.write(h1)
file.write(h2)
file.write(h3)
file.write(h4)
file.write(h5)
file.close()
c1.export(name+".svg"); print name+".svg created"
# h1.toFile(name+".txt"); print name+".txt created"
# sys.exit(0)