pythia6_qcd_jets.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 FileMC   # import FileMC
from java.util import ArrayList 
from hep.physics.vec import BasicHepLorentzVector 
from jhplot.utils import FileList
from java.awt import Color,Font
from jhplot import  HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from hephysics.particle import LParticle 
from hep.physics.jet import DurhamJetFinder
from hepsim import HepSim
import sys

url="";  Nfiles=10 # 10 files to process for checking 
default_www="http://mc.hep.anl.gov/asc/hepsim/events/ee/250gev/pythia6_qcd_all/"
if len(sys.argv)>1:
   if sys.argv[1].startswith("http"):
              flist=HepSim.getList(sys.argv[1])
              url=sys.argv[1]
   else:     # from local directory 
       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 "Found "+str(len(flist))+" files. Nr files to process= ",Nfiles

# print "Print analyse list=",flist

fjet=DurhamJetFinder(0.01) # DurhamJet with y_cut=0.01 

h1= H1D("PT(jet)",40,0,200.)  
h2= H1D("Rapidity(jet)",40,-5,5.)
#h1.doc()               # check its methods

cross=0; nev=0;  # Nfiles=len(flist)
for m in range(Nfiles):             # loop over all files in this list    
   file=FileMC(url+flist[m])        # get input file
   header = file.getHeader()
   un=float(header.getMomentumUnit()) # conversion units
   lunit=float(header.getLengthUnit())
   if m==0:
       print "ProMC v=",file.getVersion(), "E(varint)=",un,"L(varint)t=",lunit 
   for i in range(file.size()):
      nev=nev+1
      if (nev%500==0): print "Event=",nev 
      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 
      fjet.setEvent(particles)
      alljets=[] # make a new list with jets 
      for j in range(fjet.njets()):
         # print "Jet=",i," Nr of particles in jet =",fjet.nParticlesPerJet(i) 
         pjet=fjet.jet(j)   # make  HepLorentzVector
         ee=pjet.t()  # magnitude
         p3=pjet.v3() # 3-vector  
         pl=LParticle(p3.x(),p3.y(),p3.z(),ee) # convert to a class with convinient kenematics 
         alljets.append(pl)  
      for j in range(len(alljets)):
          jet=alljets[j]
          h1.fill(jet.perp())
          h2.fill(jet.rapidity())
 
   stat = file.getStat()
   cross=stat.getCrossSectionAccumulated()
   erro=stat.getCrossSectionErrorAccumulated();
   file.close()

lumi=nev/cross;
print "Lumi=%.3e pb"%lumi
print "Total cross section (pb)=",cross
c1 = HPlot("=HepSim=",600,600,1,2)
c1.visible(True)
# c1.setAutoRange()
# c1.setLogScale(1,True)
c1.setMarginLeft(100)
c1.setNameX("p_{T} [GeV]")
c1.setNameY("d σ / d p_{T} [pb/GeV]")
c1.setRange(0,200,0.0,7)

Z=h1.getDividedByBinWidth()
Z.scale(1.0/lumi)
xsec1=P1D(Z)
xsec1.setErr(1)       # show errors
xsec1.setColor(Color.blue)
c1.draw(Z)
# c1.draw(xsec1)
# h1.toTable() 

l1=HLabel("Durham jets  σ=%.3e pb"%cross, 0.3, 0.75, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l1)
l2=HLabel("=HepSim=",0.4,0.87, "NDC")
l2.setColor(Color.gray)
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l2)

# plot rapidity
c1.cd(1,2)
c1.setNameX("Rapidity")
c1.setNameY("Events")
c1.setMarginLeft(100)
c1.setAutoRange()
c1.setRangeX(-5,5)
c1.draw(h2)

# create file/image using name of the file
name="output"
if len(sys.argv[0])>0: name=sys.argv[0].replace(".py","")
file=HBook(name+".jdat","w"); print name+".jdat created"
file.write(xsec1)
file.close()
c1.export(name+".svg");    print name+".svg created"
# xsec.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)