gev380_eetunes_durham_jets.py (raw text file)

# Validation script of =HepMC= : https://atlaswww.hep.anl.gov/hepsim/
# Run it inside hs-tools or Jas4pp framework  
# 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 hep.physics.jet import DurhamJetFinder
from hep.physics.vec import BasicHepLorentzVector 
from hephysics.particle import LParticle 
import sys

fjet=DurhamJetFinder(0.02) # DurhamJet with y_cut=0.02 

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

# default data location 
default_www="http://mc.hep.anl.gov/asc/hepsim/events/ee/380gev/pythia8_tunes/"

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

# extract files with in certain tune ranges
xmin=HepSim.getRanges(flist,"_tune","_");

# reduce range
# del xmin[-1]; del xmin[-1]; 
print "Tune ranges: ",xmin

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


print "Total number of e+e- tunes =",len(filelist)
Total=len(filelist)
for ipt in range(Total):
  ptlist=filelist[ipt]
  print ipt, "Tune = "+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=",ptlist[m]
    header = file.getHeader()
    un=float(header.getMomentumUnit()) # conversion units
    for i in range(file.size()):
      nev=nev+1
      if (nev%200==0):
           if (m==0): print "Event=",nev
           else: print "Event=",nev," done=",int((100.0*m)/Nfiles),"%"
      eve = file.read(i)
      pa = eve.getParticles()    # particle information
      ve = eve.getEvent()        # event information
      particles=ArrayList()      # list of particles 
      for j in range(pa.getPxCount()): 
        if (pa.getStatus(j)==1):     # stable
           px=pa.getPx(j)/un         # pX
           py=pa.getPy(j)/un         # pY
           pz=pa.getPz(j)/un         # pZ
           e= pa.getEnergy(j)/un     # energy
           mass= pa.getMass(j)/un    # mass
           p=BasicHepLorentzVector(e,px,py,pz)
           particles.add(p) # add particle to the list 
      fjet.setEvent(particles)
      alljets=[] # make a new list with jets 
      for j in range(fjet.njets()):
         # print "Jet=",i," Nr of particles in jet =",fjet.nParticlesPerJet(i) 
         pjet=fjet.jet(j)   # make  HepLorentzVector
         ee=pjet.t()  # magnitude
         p3=pjet.v3() # 3-vector  
         pl=LParticle(p3.x(),p3.y(),p3.z(),ee) # convert to a class with convinient kenematics 
         if (pl.perp()>50): alljets.append(pl)  
      for j in range(len(alljets)):
          jet=alljets[j]
          histo[ipt].fill(jet.perp())
    stat = file.getStat()
    cross=cross+stat.getCrossSectionAccumulated() 
    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("p_{T} [GeV]")
c1.setNameY("d σ / d p_{T} [pb / GeV]")
c1.visible(True)
c1.setMarginLeft(100)
c1.setAutoRange()
c1.setRangeX(0,200)
# c1.setRangeY(0,0.1)
c1.draw(histo)
# h1.toTable() 

l1=HLabel("Pythia8 QCD σ=%.3e pb"%crossD, 0.5, 0.75, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l1)
l2=HLabel("=HepSim=",0.8,0.84, "NDC")
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
l2.setColor(Color.gray)
c1.add(l2)

l2=HLabel("Durham y=0.2",0.5,0.64, "NDC")
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)