gamma_jetphox.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
from hepsim import HepSim
import os,sys
url=""; TotalEvents=1000000000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/gamma_jetphox/"
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
nPDF=52
histoD=[] # histograms with different PDF for direc
histoF=[] # histograms with different PDF for fragm
bins=[1000,1100,1200,1300,1400,1500,1600,1800,2000,2200,2400,2600,3000,3500,4000]
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(flist)
for m in range(Nfiles): # loop over all files in this list
# if (m>5): break; # for testing only
file=FileNLO(url+flist[m]) # get input file
fragm=False
if flist[m].find("ggo_")>-1:
fragm=True
nfilesF=nfilesF+1
else: nfilesD=nfilesD+1
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:
print "ProMC v=",file.getVersion(), "M unit=",un,"L unit=",lunit
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%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)
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>1000 and Math.abs(eta)<2.5):
for k in range(nPDF):
delta = wi*(1.0-(ve.getIdata(k)/1000.)) # PDF errors are codded as (1-PDF(i)/PDF(0))*1000
if (fragm): histoF[k].fill(p.perp(),norm*delta)
else: histoD[k].fill(p.perp(),norm*delta)
file.close()
crossD=histoD[0].sumAllBinHeights()/nfilesD
#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):
newD=histoD[k].getDividedByBinWidth()
newD.scale(1.0/nfilesD)
newD.scaleErrors(0) # ignore stat. errors for weighted histograms
#newF=histoF[k].getDividedByBinWidth()
#newF.scale(1.0/nfilesF)
#newF.scaleErrors(0)
#hh=newD.oper(newF,"+")
data.append(P1D(newD))
xsec=P1D("Jetphox CT10")
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
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) # show errors
xsec.setColor(Color.blue)
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=1000
xmax=4000
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.00000011,1)
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}(γ) [GeV]")
# c1.setNameXpos(0.55, 0.65)
c1.setNameXpos(0.55, 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)
# 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)