HepSim

Repository with Monte Carlo simulations for particle physics

gamma_pythia8_ptbins.py (raw text file)
# Part of =HepSim=: http://atlaswww.hep.anl.gov/hepsim/
# Jython validation script using DataMelt http://jwork.org/dmelt/ 
# S.Chekanov (ANL)
from java.lang import *
from proto import FileMC   # import FileMC
from java.awt import Color,Font
from jhplot import  HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
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 KTjet
import math,sys

NMax=2 # max number of files per pT. Set to -1 to run over all files
default_www="http://faxbox.usatlas.org/group/hepsim/events/pp/13tev/gamma_pythia8_ptbins/"


# define Pi
Pi=math.pi
Pi2=2*math.pi

# isolation cut as in Atlas
def getIsoCut(pt):
    return 4.8+4.2E-03*pt;


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

# extract files with in certain pT ranges
"""
xmin=set()
xmax=set()
for j in flist:
      d1=j.split("_pt"); d2=d1[1].split("_") 
      xmin.add(int(d2[0]))
      xmax.add(int(d2[1]))
xmin=list(xmin); xmin.sort() 
xmax=list(xmax); xmax.sort()
"""
xmin=HepSim.getRanges(flist,"_pt","_");

# define histogram bins in ranges of pT, with 10% shift to avaoid biases
bins=[]
for j in xmin:
       bins.append(j+0.1*j)
print "Histogram bins:",bins 
  
filelist=[] # create file list using pT ranges 
histo=[]
for i in xmin:
     pt=[]
     for j in flist:
        if (j.find("pt"+str(i)+"_")>-1): pt.append(j)
     filelist.append(pt)
     histo.append(  H1D("Gamma pt"+str(i),bins) ) 

print "Total number of pT ranges=",len(filelist)
Total=len(filelist)
for ipt in range(Total):
  ptlist=filelist[ipt]
  print ipt, ") pT> "+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
    header = file.getHeader()
    un=float(header.getMomentumUnit()) # conversion units
    lunit=float(header.getLengthUnit())
    if ipt==0 and m==0:
       print "ProMC v=",file.getVersion(), "M unit=",un,"L unit=",lunit
    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

      gammas=[]
      for j in range(pa.getPxCount()):
         if (pa.getPdgId(j)==22):
           p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
           gammas.append(p)

      for p in gammas:
       esum=0
       pt_gam=p.perp()
       eta_gam=p.pseudoRapidity()
       abseta_gam=Math.abs(eta_gam)
       phi_gam=p.phi()
       for j in range(pa.getPxCount()):
          if (pa.getPdgId(j)!=22):
           g=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
           t_pt=g.perp()
           dphi=phi_gam-g.phi()
           deta=eta_gam-g.pseudoRapidity()
           adphi=Math.abs(dphi);
           if (adphi>Pi): adphi=Pi2-adphi
           ddr=Math.sqrt(deta*deta+adphi*adphi);
           # print ddr 
           if (ddr0: print "Isolation=",esum 
       if (esum

HEP.ANL.GOV