pythia8_zbbar_mbins_cross.py (raw text file)

# Validation script of =HepMC= : 
# https://atlaswww.hep.anl.gov/asc/hepsim/
# Run it inside DMelt 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
import sys

# calculate jet-jet mass 
ktjet=JetN2(0.4,"antikt",50) # anti-KT,E-mode, R=0.4,min pT=50 GeV
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.visible(True)
c1.setMarginLeft(100)
c1.setAutoRange()
c1.setRangeX(0,6000)
c1.setRangeY(0.000000001,1)
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)