HepSim

Repository with Monte Carlo simulations for particle physics

gamma_jetphox.py (raw text file)
# HepSim script using DataMelt http://jwork.org/dmelt
# Part of =HepMC= : https://atlaswww.hep.anl.gov/asc/hepsim/
# S.Chekanov (ANL)
from java.lang import *
from proto import FileNLO  # import FileNLO
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
import os,sys

url="";  TotalEvents=1000000000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/gamma_jetphox/"
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


nPDF=52
histoD=[] # histograms with different PDF for direc 
histoF=[] # histograms with different PDF for fragm 
bins=[1000,1100,1200,1300,1400,1500,1600,1800,2000,2200,2400,2600,3000,3500,4000]
for i in range(nPDF):
      h1=H1D("PT(direct) pdf="+str(i),bins)
      histoD.append(h1)
      h2=H1D("PT(fragm) pdf="+str(i),bins)
      histoF.append(h2)

nfilesF=0; nfilesD=0; 
nev=0;  Nfiles=len(flist)
for m in range(Nfiles):              # loop over all files in this list    
#   if (m>5): break;                 # for testing only 
   file=FileNLO(url+flist[m])            # get input file
   fragm=False
   if flist[m].find("ggo_")>-1:
                    fragm=True
                    nfilesF=nfilesF+1
   else: nfilesD=nfilesD+1 

   header = file.getHeader()
   un=float(header.getMomentumUnit()) # GeV units
   lunit=float(header.getLengthUnit())
   stat = file.getStat()
   cross=stat.getCrossSectionAccumulated()
   nevents=stat.getNTried()
   lumi=nevents/cross
   norm=1.0/lumi
   if m==0:
      print "ProMC v=",file.getVersion(), "M unit=",un,"L unit=",lunit

   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%20000==0): 
           if (Nfiles==1): 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
      wi=ve.getFdata(0)
      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)
           pt=p.perp()
           eta=p.pseudoRapidity()
           if (pt>1000 and Math.abs(eta)0): up=up+d*d
        if (d0: name=sys.argv[0].replace(".py","")

## save in a few formats 
file=HBook(name+".jdat","w"); print name+".jdat created"
file.write(xsec)
file.write(errors)
file.close()
c1.export(name+".svg");    print name+".svg created"
# xsec.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)

HEP.ANL.GOV