qcd_herwigpp_pt2700.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 java.awt import Color,Font
from jhplot import  HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from proto import FileMC   # import FileMC
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_herwigpp_full/qcd_herwigpp_pt2700/"
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)",50,2000,7000) 
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()

print "Start calculations.."
cross=0; nev=0;  Nfiles=len(flist)
for m in range(Nfiles):              # loop over all files in this list    
   print url+flist[m]
   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%100==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):
                 h1.fill(jets[0].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}(lead. jet) [GeV]")
c1.setNameY("d σ / d p_{T} [pb / GeV]")
c1.visible(True)
c1.setMarginLeft(90)
c1.setAutoRange(True)
#c1.setRangeX(0,10000)
#c1.setRangeY(0,10)
c1.draw(h1)
# h1.toTable() 

l1=HLabel("Herwig++ σ=%.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","")
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)