pythia8_wz_mjj.py (raw text file)

# HepSim script using DataMelt http://jwork.org/dmelt
# Part of =HepMC= : https://atlaswww.hep.anl.gov/hepsim/
# S.Chekanov (ANL)

"""
Task: Build an invariant mass of 2 jets with pT>500 GeV
      W and Z are boosted.
"""

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 hephysics.particle import LParticle 
from hepsim import HepSim
from java.util import ArrayList
from hephysics.hepsim import PromcUtil
from hephysics.jet import SCJet
import sys

url="";  TotalEvents=100000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/8tev/wzdouble_pythia/"
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("E(j)",50,500,2500) 
h2= H1D("M(jj)",50,500,3000)
#h1.setFill(True)
#h1.doc()               # check its methods

ktjet=SCJet(0.8,1,0,500,False) # R=0.6, E-mode, C-A jets,min pT=500 GeV, fast mode 
ktjet.setDebug(False)

cross=0; nev=0;  Nfiles=len(flist)
weightSumTot=0
across=0
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
   weightSum=0
   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," done=",int((100.0*m)/Nfiles),"%"
      eve = file.read(i)
      event = eve.getEvent()       # get event information  
      pa = eve.getParticles()      # particle information
      weight=  event.getWeight()   # event weight info.weight()
      weightSum = event.getPDF1()  # ad-hoc to weightSum() 
      mergingWeight = event.getPDF2() # ad-hoc to store mergingWeight()
      ptHat=event.getScale()          # pythia.info.pTHat()
      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.e())
      if (len(jets)>1):
            if (abs(jets[0].phi()-jets[1].phi())>1.5): 
                if (jets[0].perp()>500 and jets[1].perp()>500):
                  if (jets[0].mass()>50 and jets[0].mass()<120):
                    if (jets[1].mass()>50 and jets[1].mass()<120):
                         jets[0].add(jets[1])
                         h2.fill( jets[0].mass() )

   stat = file.getStat()
   across=across+stat.getCrossSectionAccumulated()
   weightSumTot=weightSumTot+weightSum
   mm=mm+1
   erro=stat.getCrossSectionErrorAccumulated();
   file.close()

cross=across/mm # average cross section 
lumi=nev/cross;
print "Lumi=%.3e pb"%lumi
print "Total cross section (pb)=",cross
print "Total weight=",weightSumTot
h1.scale(cross/(h1.getBinSize()*weightSumTot))
h2.scale(cross/(h1.getBinSize()*weightSumTot))

c1 = HPlot("HepSim",600,600,1,2) # plot histogram
c1.setNameX("E_{ j} [GeV]")
c1.setNameY("d σ / d E [pb / GeV]")
c1.visible(True)
c1.setMarginLeft(90)
c1.setAutoRange(True)
c1.cd(1,1)
c1.setRangeX(500,2500)
c1.setRangeY(0.0000001,0.001)
c1.setLogScale(1,1)
c1.draw(h1)

l1=HLabel("σ(WW,WZ,ZZ) =%.3e pb"%cross, 0.6, 0.75, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 12))
c1.add(l1)

l2=HLabel("=HepSim=",0.7,0.85, "NDC")
l2.setColor(Color.gray)
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l2)

# h1.toTable() 

c1.cd(1,2) 
c1.setLogScale(1,1)
c1.setNameX("M_{ jj} [GeV]")
c1.setNameY("d σ / d M_{ jj} [pb / GeV]")
c1.setMarginLeft(90) 
c1.setRangeX(1000,3000)
c1.setRangeY(0.0000001, 0.001)
h2.setColor(Color.blue)
h2.setFill(1)
c1.draw(h2)

# 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.write(h2)
file.close()
c1.export(name+".svg");    print name+".svg created"
# h1.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)