tev14_dijets_pythia8_mjj.py (raw text file)

# HepSim validation script using Jas4pp
# Task: Reconstruct invarinat mass of 2 jets
# S.Chekanov (ANL)
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

NMax=4 # max number of files for testing 
dataset="tev14pp_qcd_pythia8_weighted"
tag="" # only truth level
url=""
if len(sys.argv)>1:
   flist=FileList.get(sys.argv[1],"slcio")
   if (len(sys.argv)>2): NMax=int(sys.argv[2])
else:
   sites=HepSim.urlRedirector(dataset)
   url=sites[0]+"/"+tag
   flist=HepSim.getList(url)
if len(flist)==0: print "Error: No input file!"; sys.exit(0)
else: print "Reading "+str(len(flist))+" files. Nr of files= ",NMax

h1= H1D("PT(jet)",50,50,3000) 
h2= H1D("M(jj)",25,100,10000)
#h1.setFill(True)
#h1.doc()               # check its methods

ktjet=JetN2(0.4,"antikt",50) # E-mode, anti-KT,min pT=50 
# print ktjet.info()

cross=0; nev=0;  Nfiles=len(flist)
if (NMax>-1):  Nfiles=NMax           # limit to NMax files 
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
   weightSum=0
   for i in range(file.size()):
      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.perp(),weight)
      if (len(jets)>1):
                 pt1=jets[0].perp()
                 pt2=jets[1].perp()
                 if (abs(jets[0].phi()-jets[1].phi())>2):
                         jets[0].add(jets[1])
                         mass=jets[0].mass()
                         h2.fill(mass,weight)

   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("p_{T}(jet) [GeV]")
c1.setNameY("d σ / d p_{T} [pb / GeV]")
c1.visible(True)
c1.setMarginLeft(90)
c1.setAutoRange(True)
c1.cd(1,1)
c1.setRangeX(0,3000)
c1.setRangeY(0.00001,10000000)
c1.setLogScale(1,1)
c1.draw(h1)

l1=HLabel("Pythia8 QCD σ=%.3e pb"%cross, 0.5, 0.75, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 12))
c1.add(l1)

l2=HLabel("=HepSim=",0.8,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(100,10000)
c1.setRangeY(0.00000001,1000000)
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)