minbias_pythia8.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   
from jhplot.utils import FileList
from hepsim import HepSim
from java.util import ArrayList
from hephysics.hepsim import PromcUtil
from hephysics.jet import JetN2 
import sys

url="";  TotalEvents=60000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/minbias_a2_pythia8_l3/"
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)",20,0,200) 
h1.setFill(True)
#h1.doc()               # check its methods

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

across=0; nev=0;  Nfiles=len(flist)
mm=0;
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," files 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()
   across=across+stat.getCrossSectionAccumulated()
   erro=stat.getCrossSectionErrorAccumulated();
   file.close()
   mm=mm+1

cross=across/mm # average x-section
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(100)
c1.setAutoRange(True)
c1.setLogScale(1,1)
c1.setRangeX(0,200)
c1.setRangeY(100,100000000000.0)
c1.draw(h1)
# h1.toTable() 

l1=HLabel("MinBias A2 σ=%.3e pb"%cross, 0.55, 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)