minbias_pythia8.py (raw text file)

# Validation script for Jas4pp or DMelt
# Part of =HepSim= : https://atlaswww.hep.anl.gov/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

# get URL of this dataset using dataset name
name="tev100_pythia8_minbias_a14"
hepsim_server="http://atlaswww.hep.anl.gov/hepsim/"
import urllib2
response = urllib2.urlopen(hepsim_server+'/geturlmirrors.php?name='+name)
html = response.read()
response.close()
html=html.split(";")
url="";  TotalEvents=50000
default_www=html[0]
print "WWW:"+default_www
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)",40,0,200) 
h1.setFill(True)
#h1.doc()               # check its methods

ktjet=JetN2(0.5,"antikt",5) # anti-KT,E-mode, R=0.5,min pT=5 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}(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("Pythia8 MinBias A14 σ=%.3e pb"%cross, 0.55, 0.8, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 12))
c1.add(l1)

l2=HLabel("=HepSim=",0.8,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)