lepto6_dis.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 FileMC   # import FileMC
from jhplot.utils import FileList
from java.awt import Color,Font
from jhplot import  HPlot,HPlot3D,P1D,HLabel,H2D,H1D
from jhplot.io import HBook
from hephysics.particle import LParticle 
from hepsim import HepSim
import sys

url="";  Nfiles=2 # 2 files to process for checking 
default_www="http://mc.hep.anl.gov/asc/hepsim/events/ep/318gev/lepto6_dis/"
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 "Found "+str(len(flist))+" files. Nr files to process= ",Nfiles

#h1= H2D("Xbj-Q2",50,0.0,0.05,50,0,100,)  
#h1.doc()               # check its methods
h1= H1D("Xbj",1000,0.0001,1.0)  

cross=0; nev=0;  # Nfiles=len(flist)
for m in range(Nfiles):             # loop over all files in this list    
   file=FileMC(url+flist[m])        # get input file
   header = file.getHeader()
   un=float(header.getMomentumUnit()) # conversion units
   lunit=float(header.getLengthUnit())
   if m==0:
       print "ProMC v=",file.getVersion(), "M unit=",un,"L unit=",lunit 
   for i in range(file.size()):
      nev=nev+1
      if (nev%100==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
      W=ve.getScale()       # encoded  W
      Y=ve.getAlphaQED()    # encodded Y 
      X=ve.getAlphaQCD()    # x-Bj
      Q2=ve.getScalePDF() # Q2
      #print "Q2=",Q2," W=",W," Xbj=",X
      #h1.fill(X,Q2)
      h1.fill(X)
# Fill particles
#      for j in range(pa.getPxCount()): 
#         if (abs(pa.getPdgId(j))==11 or abs(pa.getPdgId(j))==13 or abs(pa.getPdgId(j))==15):
#             p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
#             h1.fill(p.perp())
 
   stat = file.getStat()
   cross=stat.getCrossSectionAccumulated()
   erro=stat.getCrossSectionErrorAccumulated();
   file.close()

lumi=nev/cross;
print "Lumi=%.3e pb"%lumi
print "Total cross section (pb)=",cross
# c1 = HPlot3D("=HepSim=",600,400)
c1=HPlot()
c1.setMarginLeft(100)
c1.setNameX("x_{Bj}")
c1.setNameY("Events")
c1.visible(True)
c1.setLogScale(0,True)
c1.setLogScale(1,True)
c1.setAutoRange()
c1.setRangeY(1,1000000)
# c1.setContour(1)
# c1.setAutoRange()
# c1.setLogScale(1,True)
# print h1.toTable()

# c1.draw(p1)
c1.draw(h1)
# h1.toTable() 

# create file/image using name of the file
name="output"
if len(sys.argv[0])>0: name=sys.argv[0].replace(".py","")
f=HBook(name+".jdat","w"); print name+".jdat created"
f.write(h1)
f.close()
c1.export(name+".svg");    print name+".svg created"
# xsec.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)