jets_nlojetpp.py (raw text file)# HepSim script using DataMelt http://jwork.org/dmelt
# Part of =HepMC= : https://atlaswww.hep.anl.gov/asc/hepsim/
# S.Chekanov (ANL)
from java.lang import *
from java.awt import Color,Font
from jhplot import HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from proto import FileNLO # import FileNLO
from jhplot.utils import FileList
from hephysics.particle import LParticle
import glob,os, sys
from hepsim import HepSim
from javax.swing import JOptionPane
if JOptionPane.showConfirmDialog(None, "This macro works for local files. Quit?","Warning", JOptionPane.OK_CANCEL_OPTION) == JOptionPane.OK_OPTION: System.exit(0);
url=""
if len(sys.argv)>1: # if one argument is given, read all files from a directory
DIR=sys.argv[1]
flist50=glob.glob(DIR+"*_pt50_*promc")
flist100=glob.glob(DIR+"*_pt100_*promc")
flist200=glob.glob(DIR+"*_pt200_*promc")
flist400=glob.glob(DIR+"*_pt400_*promc")
flist800=glob.glob(DIR+"*_pt800_*promc")
flist1600=glob.glob(DIR+"*_pt1600_*promc")
flist3200=glob.glob(DIR+"*_pt3200_*promc")
flist6400=glob.glob(DIR+"*_pt6400_*promc")
if (len(flist50)*len(flist100)*len(flist200)*len(flist800)*len(flist1600)*len(flist3200)*len(flist6400)==0):
print "Error: some files are missing !"; sys.exit(0)
flist= [ flist50, flist100, flist200, flist400, flist800, flist1600, flist3200, flist6400]
else:
url="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/jets_nlojetpp/"
flist=HepSim.getList(url)
if len(flist)==0: print "Error: No input file!"; sys.exit(0)
minPT=50
maxETA=3.5
print "Inc. cross section. minPT=",minPT," max Eta=",maxETA
SYST=1+7+40 # 0, central, 7- scale uncertainties, 40-PDF (MSTW)
bins=[50,100,200,400,800,1600,3200,6400,12800]
# bins=[50,75,100,150,200,300,400,600,800,1200,1600,2400,3200,6400]
nev=0
ptbins=[[50,100], [100,200],[200,400], [400,800],[800,1600],[1600,3200],[3200,6400],[6400,10000]]
pthisto=[]
print "Process all pT ranges:",ptbins
for ipt in range(len(flist)):
ptlist=flist[ipt]
print "Dealing with="+str(ptbins[ipt])+" GeV"
Nfiles=len(ptlist)
totruns=0
avcross=0
nfiles=0
histo=[] # histograms with different PDF for direct
for i in range(SYST):
h1=H1D("PT pdf="+str(i),bins)
histo.append(h1)
for m in range(len(ptlist)): # loop over all files in this list
# if (m>10): break # debugging mode
bsize=os.path.getsize(ptlist[m])
if (bsize<1000):
print ptlist[m] + " too small -> probably broken"
continue
file=FileNLO(ptlist[m]) # get input file
header = file.getHeader()
un=float(header.getMomentumUnit()) # GeV units
lunit=float(header.getLengthUnit())
stat = file.getStat()
cross=stat.getCrossSectionAccumulated()
avcross=avcross+cross
nevents=stat.getNTried()
nfiles=nfiles+1
if m==0:
print "ProMC v=",file.getVersion(), "M unit=",un," L unit=",lunit
for i in range(file.size()):
nev=nev+1
if (nev%20000==0):
if (Nfiles==1): print "Event=",nev
else: print "Event=",nev," done=",int((100.0*m)/Nfiles),"%"
eve = file.read(i)
pa = eve.getParticles() # particle information
ve = eve.getEvent() # event information
wi=ve.getFdata(0) # event weight
runs=ve.getProcessID() # runs (not events!)
jets=[]
for j in range(pa.getPxCount()):
p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
pt=p.perp()
eta=p.pseudoRapidity()
if (pt>ptbins[ipt][0] and pt0): up=up+d*d
if (d<0): dw=dw+d*d
xsec.add(x,y,errX,errX,Math.sqrt(up),Math.sqrt(dw))
xsec.setErr(0) # make sure you show errors
xsec.setErrSys(0)
pthisto.append(xsec)
# add all pT bins
xsec=pthisto[0]
xsec.setTitle("Cross section")
for i in range(1,len(flist)): xsec.operPlusY(pthisto[i])
errors=P1D("Relative errors")
for j in range(xsec.size()):
y=xsec.getY(j)
if y !=0:
errors.add(xsec.getX(j),1,xsec.getXleft(j),xsec.getXright(j),xsec.getYupper(j)/y,xsec.getYlower(j)/y)
errors.setErr(1)
# errors.setErrFillColor(Color.blue,0.3)
c1 =HPlot("=HepSim=",400,600, 1,2)
c1.visible(True)
c1.resizePad(0,0,1.0, 1.5)
c1.resizePad(0,1,1.0, 0.5)
xmin=50
xmax=12800
c1.cd(1,1)
c1.setNameY("d σ / d p_{T} [pb/GeV]")
c1.setNameYpos(-0.01, 0.93)
c1.setMarginLeft(55)
c1.setMarginBottom(0.0)
c1.setTicLabels(0, False)
c1.setRange(xmin,xmax,0.000000011,1000000)
c1.setLogScale(1,1)
tot=0
for i in range(xsec.size()):
tot=tot+xsec.getY(i)*(bins[i+1]-bins[i])
l1=HLabel("NLO σ=%.3e pb"%tot, 0.3, 0.7, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 15))
c1.add(l1)
l2=HLabel("=HepSim=",0.5,0.85, "NDC")
l2.setColor(Color.gray)
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l2)
c1.draw(xsec)
c1.cd(1,2)
c1.setRange(xmin,xmax,0.7,1.3)
c1.setLegend(False)
c1.setNameY("Ratio")
c1.setNameYpos(-0.01, 0.9)
c1.setMarginTop(0.0)
c1.setMarginLeft(55)
c1.setNameX("p_{T}(jet) [GeV]")
c1.setNameXpos(0.55, 0.65)
from jhplot.shapes import Line
line = Line(xmin,1.0, xmax, 1.0)
line.setPosCoord("USER")
line.setTransparency(0.5)
c1.add(line)
c1.draw(errors)
# create file/image using name of the file
name="output"
if len(sys.argv[0])>0: name=sys.argv[0].replace(".py","")
## save in a few formats
file=HBook(name+".jdat","w"); print name+".jdat created"
file.write(xsec)
file.write(errors)
file.close()
c1.export(name+".svg"); print name+".svg created"
# xsec.toFile(name+".txt"); print name+".txt created"
sys.exit(0)