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 JetN2
import sys
url=""; TotalEvents=20000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/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("M(j)",100,0,500)
h2= H1D("M(jj)",80,2000,20000)
#h1.setFill(True)
#h1.doc() # check its methods
ktjet=JetN2(0.8,"antikt",500) # anti-KT,E-mode, R=0.8,min pT=500 GeV
print ktjet.info()
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.mass())
if (len(jets)>1):
if (abs(jets[0].phi()-jets[1].phi())>2):
if (jets[0].perp()>1000 and jets[1].perp()>1000):
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("M_{ j} [GeV]")
c1.setNameY("d σ / d M [pb / GeV]")
c1.visible(True)
c1.setMarginLeft(90)
c1.setAutoRange(True)
c1.cd(1,1)
c1.setRangeX(0,500)
c1.setRangeY(0.00001,0.1)
c1.setLogScale(1,1)
c1.draw(h1)
l1=HLabel("σ(WW,WZ,ZZ) =%.3e pb"%cross, 0.5, 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(2000,20000)
c1.setRangeY(0.00000001, 0.1)
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)