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.awt import Color,Font
from jhplot import  HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from proto import FileMC   # 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
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/100tev/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=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("Higgs 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 range = "+str(xmin[ipt])+"-"+str(xmax[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