HepSim

Repository with Monte Carlo simulations for particle physics

ttbar_pythia8_ptbins.py (raw text file)
# Part of =HepSim=: http://atlaswww.hep.anl.gov/hepsim/
# Jython validation script using DataMelt http://jwork.org/dmelt/ 
# S.Chekanov (ANL)

from java.lang import *
from proto import FileMC   # import FileMC
from java.awt import Color,Font
from jhplot import  HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from jhplot.utils import FileList
from hephysics.particle import LParticle 
from hepsim import HepSim
from java.util import ArrayList
from hephysics.hepsim import PromcUtil
import sys

NMax=5 # max number of files per pT. Set to -1 to run over all files
# default data location 
default_www="http://faxbox.usatlas.org/group/hepsim/events/pp/100tev/ttbar_pythia8_ptbins/"

# get data either from local file system or stream them from URL
url="" 
TotalEvents=100000 # number of events per sample 
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

# extract files with in certain pT ranges
xmin=HepSim.getRanges(flist,"_pt","_");

# define histogram bins in ranges of pT, with 10% shift to avaoid biases
bins=[]
for j in xmin:
       bins.append(j)
print "Histogram bins",bins 

# create file list using pT ranges 
filelist=[]
histo=[]
for i in xmin:
     pt=[]
     for j in flist:
        if (j.find("pt"+str(i)+"_")>-1): pt.append(j)
     filelist.append(pt)
     histo.append(  H1D("pt"+str(i),bins) ) # 50 GeV bins 

print "Total number of pT ranges=",len(filelist)
Total=len(filelist)
for ipt in range(Total):
  ptlist=filelist[ipt]
  print ipt, ") Min pt = "+str(xmin[ipt]) 
  cross=0; xcross=0; nev=0;  Nfiles=len(ptlist)
  if (NMax>-1):  Nfiles=NMax           # limit to NMax files 
  for m in range(Nfiles):              # loop over all files in this list    
    file=FileMC(url+ptlist[m])         # get input file
    header = file.getHeader()
    un=float(header.getMomentumUnit()) # conversion units
    lunit=float(header.getLengthUnit())
    if ipt==0 and m==0:
       print "ProMC v=",file.getVersion(), "M unit=",un,"L unit=",lunit
    print "Read=",url+ptlist[m]
    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
      nev=nev+1
      if (nev%200==0):
           if (Nfiles==1): print "Event=",nev
           else: print "Event=",nev," done=",int((100.0*m)/Nfiles),"%", 'x-cross=',xcross 
      eve = file.read(i)
      pa = eve.getParticles()    # particle information
      ve = eve.getEvent()        # event information
      top=[];  atop=[]
      for j in range(pa.getPxCount()):
           p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un,pa.getMass(j)/un)
           if (pa.getPdgId(j)==6):   top.append(p)
           if (pa.getPdgId(j)==-6): atop.append(p)
           if len(top)==1 and len(atop)==1:  break  # take only first
      if (len(top)==1 and len(atop)==1):
              histo[ipt].fill(top[0].perp())
              histo[ipt].fill(atop[0].perp())
    stat = file.getStat()
    xcross=stat.getCrossSectionAccumulated()
    cross=cross+xcross;
    erro=stat.getCrossSectionErrorAccumulated();
    file.close()
  cross=cross/Nfiles
  lumi=nev/cross
  histo[ipt].scale(1.0/lumi)   # for cross section 
  print ipt, ") Cross section = %.3e pb"%cross, " Lumi=%.3e pb"%lumi

h1=histo[0] #  combine all pT ranges 
h1.setTitle("t,tbar cross section")
for i in range(1,len(histo)):
    h1=h1.oper(histo[i],"+")

crossD=h1.sumAllBinHeights() # get total cross section 
h1=h1.getDividedByBinWidth() # get differential cross section 
c1 = HPlot("HepSim",600,400) # plot histogram
c1.setNameX("p_{T}(t) [GeV]")
c1.setNameY("d σ / d p_{T} (t) [pb / GeV]")
c1.visible(True)
c1.setMarginLeft(100)
c1.setAutoRange()
c1.setRangeX(0,30000)
c1.setRangeY(10e-16,10e+6)
c1.setLogScale(1,True)
c1.draw(h1)
# h1.toTable() 

l1=HLabel("Pythia8   QCD σ=%.3e pb"%crossD, 0.5, 0.65, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l1)
l2=HLabel("=HepSim=",0.65,0.84, "NDC")
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
l2.setColor(Color.gray)
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)

HEP.ANL.GOV