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 proto import FileNLO # import FileNLO
from java.awt import Color,Font
from jhplot import HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from jhplot.utils import FileList
from hephysics.particle import LParticle
import glob,os, sys
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=None
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")
print len(flist50),len(flist100),len(flist200),len(flist800),len(flist1600),len(flist3200)
if (len(flist50)*len(flist100)*len(flist200)*len(flist800)*len(flist1600)*len(flist3200)==0):
print "Error: some files are missing !"; sys.exit(0)
else:
url="http://atlaswww1.hep.anl.gov/asc/RefHepSim/events"
flist=[url+"/pp/8tev/jets_nlojetpp/nlojetpp_001.promc"]
print "Reading 1 file from WWW=",flist[0]
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=7+40 # 0, central, 7- scale uncertainties, 40-PDF (MSTW)
bins=[50,100,200,400,800,1600,3200,6400]
# bins=[50,100,200,300,400,600,800,1200,1600,3200,6400]
# bins=[50,75,100,150,200,300,400,600,800,1200,1600,2400,3200,6400]
# bins=[50,75,100,150,200,300,400,600,800,1000,1200,1400,1600,2000,2400,3200,4000]
nev=0
ptbins=[[50,100], [100,200],[200,400], [400,800],[800,1600],[1600,3200],[3200,6400]]
flist= [ flist50, flist100, flist200, flist400, flist800, flist1600, flist3200]
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):
histo.append( H1D("PT ="+str(i),bins) )
for m in range(Nfiles): # loop over all files in this list
# if (m>1): 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
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): scale_up=scale_up+d*d
if (n<7 and d<0): scale_dw=scale_dw+d*d
if (d>0): up=up+d*d
if (d<0): dw=dw+d*d
xsec.add(x,y,errX,errX,Math.sqrt(scale_up),Math.sqrt(scale_dw),0,0,Math.sqrt(up),Math.sqrt(dw))
# print xsec.toString()
pthisto.append(xsec)
# add all Y-values bins
pp=pthisto[0]
for i in range(1,len(flist)): pp.operPlusY(pthisto[i])
# get data with errors
errors=P1D("Relative errors")
errors.setErr(1) # make sure you show errors
errors.setErrSys(1)
for j in range(pp.size()):
y=pp.getY(j)
if y !=0:
errors.add(pp.getX(j),1,pp.getXleft(j),pp.getXright(j),pp.getYupper(j)/y,pp.getYlower(j)/y,0,0,pp.getYupperSys(j)/y,pp.getYlowerSys(j)/y)
pp.setTitle("PT(γ)")
pp.setColor(Color.blue)
pp.setErr(1) # make sure you show errors
pp.setErrSys(1)
pp.setTitle("PT(jet)")
tot=0
for i in range(pp.size()):
tot=tot+pp.getY(i)*(bins[i+1]-bins[i])
print "Cross section=",tot," pb"
c1 =HPlotJa("Canvas",400,600, 1,2)
c1.setAntiAlias(False)
c1.visible()
xmin=0
xmax=6000
c1.cd(1,1)
c1.showKey(False)
c1.setLabelShift(1, 0.0, 0.27)
lab=c1.getLabel(1)
lab.setText("d σ / d p_{T}(jet) [pb / GeV]")
lab=c1.getLabel(0)
lab.setText("")
c1.setShowStatBox(False)
c1.setTicksLabels(0, False)
c1.showKey(False)
c1.setPadSize(0.79,0.68)
c1.setPadLocation(0.16,0.04)
c1.setRange(xmin,xmax,0.00000001,10E05)
c1.setLogScale(1,1)
#p1.setErrFillColor(Color.yellow,0.3)
c1.draw(pp)
c1.cd(1,2)
c1.showKey(False)
c1.setLabelShift(1, 0.0, 0.22)
labY=c1.getLabel(1)
labY.setText("Ratio")
c1.setLabelShift(0, -0.2, 0)
labX=c1.getLabel(0)
labX.setText("p_{T}(jet) [GeV]")
f=F1D("1",xmin,xmax)
f.setPenWidth(1)
c1.draw(f)
c1.setShowStatBox(False)
c1.setRange(xmin,xmax,0.4,1.6)
c1.setPadSize(0.79,0.18)
c1.setPadLocation(0.16, 0.73)
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","")
c1.export(name+".png"); print name+".png created"
xsec.toFile(name+".txt"); print name+".txt created"
## save for future
file=HBook(name+".jdat","w"); print name+".jdat created"
file.write(pp)
file.write(errors)
for k in range(len(flist)):
file.write(data[k])
file.close()