qcd_pythia8_pt8000.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.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 hephysics.particle import LParticle 
from hepsim import HepSim
from java.util import ArrayList
from hephysics.hepsim import PromcUtil
from hephysics.jet import JetN2 
import sys

url="";  TotalEvents=10000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/qcd_pythia8_pt/qcd_pythia8_pt8000/"
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

h1= H1D("PT(jet lead)",200,5000,20000) 
h1.setFill(True)
#h1.doc()               # check its methods

ktjet=JetN2(0.5,"antikt",100) # anti-KT,E-mode, R=0.5,min pT=100 GeV,fast
print ktjet.info()

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
   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 # stop after some events
      nev=nev+1
      if (nev%200==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
      particles=PromcUtil.getParticleDList(file.getHeader(), pa, 1, -1, 1000);
      ktjet.buildJets(particles);
      # ktjet.printJets();
      jets=ktjet.getJetsSorted()
      if (len(jets)>0):
                 lead=jets[0]
                 h1.fill(lead.perp())
   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
h1.scale(1.0/(h1.getBinSize()*lumi))
c1 = HPlot("HepSim",600,400) # plot histogram
c1.setNameX("p_{T}(jet) [GeV]")
c1.setNameY("d σ / d p_{T} [pb / GeV]")
c1.visible(True)
c1.setMarginLeft(100)
c1.setAutoRange(True)
#c1.setRangeX(0,10000)
#c1.setRangeY(0,10)
c1.draw(h1)
# h1.toTable() 

l1=HLabel("Pythia8 QCD σ=%.3e pb"%cross, 0.5, 0.8, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 12))
c1.add(l1)

l2=HLabel("=HepSim=",0.7,0.87, "NDC")
l2.setColor(Color.gray)
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l2)

# 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"
file.write(h1)
file.close()
c1.export(name+".svg");    print name+".svg created"
# h1.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)