pythia8_zbbar_mbins.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.awt import Color,Font
from java.lang import *
from proto import FileMC
from jhplot import  HPlot,P1D,HLabel,H1D 
from jhplot.utils import FileList
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.4,"antikt",50) # anti-KT,E-mode, R=0.5,min pT=100 GeV,fast
print ktjet.info()

NMax=2 # max number of files per mass. Set to -1 to run over all files

# default location 
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/13tev/pythia8_zbbar_nnpdf/"

# get data either from local file system or stream them from URL
url="" 
TotalEvents=10000 # 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

xmin=HepSim.getRanges(flist,"_m","_");
print "Mass ranges: ",xmin

# create file list using mass ranges 
filelist=[]
histo=[]
n=0
for i in xmin:
     ma=[]
     for j in flist:
        if (j.find("m"+str(i)+"_")>-1): ma.append(j)
     filelist.append(ma)
     hh=H1D("M="+str(i),100,0,6000)
     hh.setColor(Color(50+n*20,20+n*10,200-5*n)) # set color map 
     histo.append( hh ) # 50 GeV bins 
     n=n+1

print "Total number of mass ranges=",len(filelist)
Total=len(filelist)
for ipt in range(Total):
  ptlist=filelist[ipt]
  print ipt, "Mass = "+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()
      # mass of 2 leading jets
      if (len(jets)>1):
                 ptlead1=jets[0].perp()
                 ptlead2=jets[1].perp()
                 #print "jet1=",jets[0]
                 #print "jet2=",jets[1]
                 if (abs(jets[0].phi()-jets[1].phi())>2):
                         jets[0].add(jets[1])
                         mass=jets[0].mass()
                         # print "mass=",mass
                         histo[ipt].fill(mass)

    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 
  histo[ipt]=histo[ipt].getDividedByBinWidth() # get differential cross section 
  # print ipt, ") Cross section = %.3e pb"%cross, " Lumi=%.3e pb"%lumi

crossD=histo[0].sumAllBinHeights() # get total cross section 
c1 = HPlot("HepSim",600,400) # plot histogram
c1.setNameX("M_{ jj} [GeV]")
# c1.setNameY("d σ / d M_{jj} [pb / GeV]")
c1.setNameY("Events / GeV")
c1.visible(True)
c1.setMarginLeft(100)
c1.setAutoRange()
c1.setRangeX(0,6000)
c1.setRangeY(0,30)
# c1.setLogScale(1,True)
# plot for mass=4 TeV
c1.draw(histo)
# 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.75,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","")
from jhplot.io import HBook
file=HBook(name+".jdat","w"); print name+".jdat created"
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)