gamma_jetphox_ptbins.py (raw text file)# HepSim script using DataMelt http://jwork.org/dmelt
# Shows gamma cross section with 52 PDF uncertainties
# S.Chekanov (ANL)
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 FileNLO
from jhplot.utils import FileList
from hephysics.particle import LParticle
import os, sys
import glob
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);
flist=[] # file list
ptbins=[[200,1500]]
url=None
print "Nr of arguments=",len(sys.argv)
if len(sys.argv)>1: # if one argument is given, read all files from a directory
DIR=sys.argv[1]
print DIR
flist200=glob.glob(DIR+"*_mu1_200_600_*promc")
flist600=glob.glob(DIR+"*_mu1_600_1200_*promc")
flist1200=glob.glob(DIR+"*_mu1_1200_2400_*promc")
flist2400=glob.glob(DIR+"*_mu1_2400_4800_*promc")
flist4800=glob.glob(DIR+"*_mu1_4800_15000_*.promc")
print len(flist200), len(flist600), len(flist1200), len(flist2400), len( flist4800)
if (len(flist200)*len(flist600)*len(flist1200)*len(flist2400)*len(flist4800)==0):
print "Error: some files are missing !"; sys.exit(0)
flist= [flist200, flist600, flist1200, flist2400, flist4800]
ptbins=[[200,600],[600,1200],[1200,2400],[2400,4800],[4800,15000]]
else:
url="http://mc.hep.anl.gov/asc/hepsim/events"
flist=[url+"/pp/100tev/gamma_jetphox/ggd_mu1_1000_6000_run0_atlas50.promc"]
print "Reading 1 file from WWW=",flist[0]
isoCone=0.4
isoET=7
EtaMax=2.5
nPDF=41
bins=[200,300,400,500,600,700,800,900,1000,1200,1400,1600,1800,2000,2200,2400,2600,3000,3400,3800,4200,4800,5000,6000,7000,8000,10000,14000]
pthisto=[]
for ipt in range(len(flist)):
ptlist=flist[ipt]
histoD=[] # histograms with different PDF for direc
histoF=[] # histograms with different PDF for fragm
print "Dealing with="+str(ptbins[ipt])+" GeV"
for i in range(nPDF):
h1=H1D("PT(direct) pdf="+str(i),bins)
histoD.append(h1)
h2=H1D("PT(fragm) pdf="+str(i),bins)
histoF.append(h2)
nfilesF=0; nfilesD=0;
nev=0; Nfiles=len(ptlist)
for m in range(len(ptlist)): # loop over all files in this list
# if (m>10): break;
bsize=os.path.getsize(ptlist[m])
if (bsize<500):
print ptlist[m] + " too small -> probably broken"
continue
file=FileNLO(ptlist[m]) # get input file
fragm=False
if ptlist[m].find("ggo_")>-1:
fragm=True
nfilesF=nfilesF+1
else: nfilesD=nfilesD+1
#print "nfilesD=",nfilesD, " nfilesF=",nfilesF
header = file.getHeader()
un=float(header.getMomentumUnit()) # GeV units
lunit=float(header.getLengthUnit())
stat = file.getStat()
cross=stat.getCrossSectionAccumulated()
nevents=stat.getNTried()
lumi=nevents/cross
norm=1.0/lumi
if m==0 and ipt==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):
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)
if (ve.getIdataCount() != nPDF): print "Error in PDF set!"
for j in range(pa.getPxCount()):
if (pa.getPdgId(j)==22):
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) :
crossD=histoD[0].sumAllBinHeights()/nfilesD
if (nfilesF>0) :
crossF=histoF[0].sumAllBinHeights()/nfilesF
print "Direct cont =",crossD," from ",nfilesD," files"
print "Fragm cont =",crossF," from ",nfilesF," files"
data=[]
for k in range(nPDF):
if (nfilesD>0) :
newD=histoD[k].getDividedByBinWidth()
newD.scale(1.0/nfilesD)
newD.scaleErrors(0)
if (nfilesF>0) :
newF=histoF[k].getDividedByBinWidth()
newF.scale(1.0/nfilesF)
newF.scaleErrors(0)
if nfilesD>0 and nfilesF>0:
hh=newD.oper(newF,"+")
data.append(P1D(hh))
if nfilesD>0 and nfilesF==0:
data.append(P1D(newD))
if (len(data)==0):
print "Error: No files!"; sys.exit(0)
# create P1D with 2nd level errors evaluated from an array of P1D
#xsec=data[0]
#xsec=xsec.getSys(data)
#xsec.setErrToZero(1) # set all stat. error to 0
xsec=P1D("Cross sections:"+str(ptbins[ipt]))
for j in range(data[0].size()):
x=data[0].getX(j)
y=data[0].getY(j)
errX=0.5*(bins[j+1]-bins[j])
up=0 # scale+PDF
dw=0
for n in range(1,nPDF):
d=data[n].getY(j)-y
if (d>0): 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(1) # make sure you show errors
xsec.setErrSys(1)
pthisto.append(xsec)
"""
print pthisto[0].toString()
print pthisto[1].toString()
print pthisto[2].toString()
print pthisto[3].toString()
"""
# merge PT bins
pp=pthisto[0]
for i in range(1,len(flist)): pp.operPlusY(pthisto[i])
pp.setTitle("NLO QCD (MSTW2008 NLO)")
# pp.setColor(Color.blue)
pp.setErr(1) # make sure you show errors
pp.setErrSys(1)
# get data with errors
errors=P1D("Relative errors")
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)
errors.setErr(1) # make sure you show errors
errors.setErrSys(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=0
xmax=15000
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.0000000011,1000)
c1.setLogScale(1,1)
l2=HLabel("=HepSim=",0.3,0.6, "NDC")
l2.setColor(Color.blue)
c1.add(l2)
c1.draw(pp)
c1.cd(1,2)
c1.setRange(xmin,xmax,0.85,1.15)
c1.setLegend(False)
c1.setNameY("Ratio")
c1.setNameYpos(-0.01, 0.9)
c1.setMarginTop(0.0)
c1.setMarginLeft(55)
c1.setNameX("p_{T}(γ) [GeV]")
# c1.setNameXpos(0.55, 0.65)
c1.setNameXpos(0.6, 0.1)
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)
tot=0
for i in range(xsec.size()):
tot=tot+xsec.getY(i)*(bins[i+1]-bins[i])
print "JETPHOX NLO |η|<2.5 σ=%.3e pb"%tot
# 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(pp)
file.write(errors)
for k in range(len(flist)): file.write(pthisto[k])
file.close()
c1.export(name+".svg"); print name+".svg created"
# xsec.toFile(name+".txt"); print name+".txt created"
# sys.exit(0)