HepSim

Repository with Monte Carlo simulations for particle physics

qcd_pythia8_ptbins_jets.py (raw text file)
# Validation script of =HepMC= : https://atlaswww.hep.anl.gov/asc/hepsim/
# Run it inside 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
from hephysics.jet import JetN2,ParticleD
import sys

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

NMax=2 # max number of files per pT. Set to -1 to run over all files
# default location 
default_www="http://faxbox.usatlas.org/group/hepsim/events/pp/13tev/qcd_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=set()
xmax=set()
for j in flist:
      d1=j.split("_pt"); d2=d1[1].split("_") 
      xmin.add(int(d2[0]))
      xmax.add(int(d2[1]))
xmin=list(xmin); xmin.sort() 
xmax=list(xmax); xmax.sort()
"""
xmin=HepSim.getRanges(flist,"_pt","_");

# reduce range
#del xmin[-1]; del xmin[-1]; 
#del xmax[-1]; del xmax[-1]; 

# 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, ") pT range = "+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
    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 (m==0): 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
      particles=PromcUtil.getParticleDList(file.getHeader(), pa, 1, -1, 1000);
      ktjet.buildJets(particles);
      # ktjet.printJets();
      jets=ktjet.getJetsSorted()
      if (len(jets)>0):
                 ptlead=jets[0].perp() 
                 if ptlead>xmin[ipt]: 
                         histo[ipt].fill(ptlead)
    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("All bins")
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}(jet lead.) [GeV]")
c1.setNameY("d σ / d p_{T} (jet lead.) [pb / GeV]")
c1.visible(True)
c1.setMarginLeft(100)
# c1.setAutoRange()
c1.setRangeX(0,7000)
c1.setRangeY(10e-15,10e+7)
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)
for i in histo: file.write(i)
file.close()
c1.export(name+".svg");    print name+".svg created"
# h1.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)

HEP.ANL.GOV