tev3_eetunes_eventshapes.py (raw text file)# Validation script of =HepMC= : https://atlaswww.hep.anl.gov/hepsim/
# Run it inside hs-tools or Jas4pp framework
# Task: Calculations of event shapes in e+e-
# S.Chekanov (ANL)
from java.awt import Color,Font
from java.lang import *
from proto import FileMC
from jhplot import HPlot,P1D,HLabel,H1D
from jhplot.utils import FileList
from hepsim import HepSim
from java.util import ArrayList
from hep.physics.vec import BasicHepLorentzVector
from hep.physics.jet import EventShape
from hephysics.particle import LParticle
import sys
evshape=EventShape()
NMax=2 # max number of files per tune. Set to -1 to run over all files
# default data location
default_www="http://mc.hep.anl.gov/asc/hepsim/events/ee/3tev/pythia8_tunes/"
# get data either from local file system or stream them from URL
url=""
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"
# extract files with in certain tune ranges
xmin=HepSim.getRanges(flist,"_tune","_");
# reduce range
# del xmin[-1]; del xmin[-1];
print "Tune ranges: ",xmin
# create file list using mass ranges
filelist=[]
histo1=[] # oblateness: Major Thrust - Minor Thrust
histo2=[] # Major Thrust
histo3=[] # Minor Thrust
n=0
for i in xmin:
ma=[]
for j in flist:
if (j.find("tune"+str(i)+"_")>-1): ma.append(j)
filelist.append(ma)
hh=H1D("tune="+str(i),40,0,0.2)
hh.setColor(Color(50+n*20,20+n*10,200-n*4)) # set color map
histo1.append( hh )
hh=H1D("tune="+str(i),40,0,0.2)
hh.setColor(Color(50+n*20,20+n*10,200-n*4)) # set color map
histo2.append( hh )
hh=H1D("tune="+str(i),40,0,0.2)
hh.setColor(Color(50+n*20,20+n*10,200-n*4)) # set color map
histo3.append( hh )
n=n+1
print "Total number of e+e- tunes =",len(filelist)
Total=len(filelist)
for ipt in range(Total):
ptlist=filelist[ipt]
print ipt, "Tune = "+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=",ptlist[m]
header = file.getHeader()
un=float(header.getMomentumUnit()) # conversion units
for i in range(file.size()):
nev=nev+1
if (nev%200==0):
if (m==0): 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
particles=ArrayList() # list of particles
for j in range(pa.getPxCount()):
if (pa.getStatus(j)==1): # stable
px=pa.getPx(j)/un # pX
py=pa.getPy(j)/un # pY
pz=pa.getPz(j)/un # pZ
e= pa.getEnergy(j)/un # energy
mass= pa.getMass(j)/un # mass
p=BasicHepLorentzVector(e,px,py,pz)
particles.add(p) # add particle to the list
# see: http://java.freehep.org/lib/freehep/api/hep/physics/jet/EventShape.html
evshape.setEvent(particles)
oblateness=evshape.oblateness()
majoraxis=evshape.majorAxis()
minoraxis=evshape.minorAxis()
trustaxis=evshape.thrustAxis()
trust=evshape.thrust()
histo1[ipt].fill(oblateness)
histo2[ipt].fill(trust.y()) # Major
histo3[ipt].fill(trust.z()) # Minor
stat = file.getStat()
cross=cross+stat.getCrossSectionAccumulated()
erro=stat.getCrossSectionErrorAccumulated()
file.close()
cross=cross/Nfiles
lumi=nev/cross
histo1[ipt].scale(1.0/lumi) # for cross section
histo1[ipt]=histo1[ipt].getDividedByBinWidth()
histo2[ipt].scale(1.0/lumi) # for cross section
histo2[ipt]=histo2[ipt].getDividedByBinWidth()
histo3[ipt].scale(1.0/lumi) # for cross section
histo3[ipt]=histo3[ipt].getDividedByBinWidth()
# print ipt, ") Cross section = %.3e pb"%cross, " Lumi=%.3e pb"%lumi
crossD=histo1[0].sumAllBinHeights() # get total cross section
c1 = HPlot("HepSim",700,800, 1, 3) # plot histogram
c1.setNameX("O")
c1.setNameY("d σ / d O [pb]")
c1.visible(True)
c1.setMarginLeft(90)
c1.setAutoRange()
c1.setRangeX(0,0.2)
# c1.setRangeY(0,0.1)
c1.draw(histo1)
# h1.toTable()
l1=HLabel("Pythia8 QCD σ=%.3e pb"%crossD, 0.5, 0.75, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l1)
l2=HLabel("=HepSim=",0.8,0.84, "NDC")
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
l2.setColor(Color.gray)
c1.add(l2)
c1.cd(1,2)
c1.setAutoRange()
c1.setMarginLeft(90)
c1.setNameX("T_{Major}")
c1.setNameY("d σ / d T_{Major} [pb]")
c1.setRangeX(0,0.2)
# c1.setRangeY(0,0.1)
c1.draw(histo2)
c1.cd(1,3)
c1.setAutoRange()
c1.setMarginLeft(90)
c1.setRangeX(0,0.2)
# c1.setRangeY(0,0.1)
c1.setNameX("T_{Minor}")
c1.setNameY("d σ / d T_{Minor} [pb]")
c1.draw(histo3)
# create file/image using name of the file
name="output"
if len(sys.argv[0])>0: name=sys.argv[0].replace(".py","")
from jhplot.io import HBook
file=HBook(name+".jdat","w"); print name+".jdat created"
for i in histo1: file.write(i)
for i in histo2: file.write(i)
for i in histo3: file.write(i)
file.close()
c1.export(name+".svg"); print name+".svg created"
# h1.toFile(name+".txt"); print name+".txt created"
# sys.exit(0)