check_mu40.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 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
import sys,math 

url="";  TotalEvents=10000000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/misc/pgun/pgun_pi_mu40/"
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

h1= H1D("dR",100,0,0.2)  
h1.setFill(True)
h2= H1D("N(charged)",30,0,30) 
h2.setFill(True)
h3= H1D("N(stable)",30,0,30)  
h3.setFill(True)

pi=3.1415
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 
   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%500==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

      # first stable particle is a gun
      eta=0
      phi=0
      for j in range(pa.getPxCount()): 
        if (pa.getStatus(j)==1):
           p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
           pmom=p.px()*p.px() + p.py()*p.py() + p.pz()*p.pz() 
           pmom=math.sqrt(pmom)
           eta=p.pseudoRapidity();
           phi=p.phi();
           break
           #print eta, phi
      nn=0
      nns=0
      # skip 1st and 2nd (gun)
      for j in range(2,pa.getPxCount()):
        if (pa.getStatus(j)==1):
           p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
           deta=p.pseudoRapidity()-eta
           dphi=p.phi()-phi
           # print pa.getCharge(j), pa.getPdgId(j)
           if (dphi>pi): dphi=2*pi-dphi
           dR=math.sqrt(deta*deta + dphi*dphi)
           h1.fill(dR)
           nns=nns+1
           if (pa.getCharge(j) !=0): nn=nn+1 
      h2.fill(nn)
      h3.fill(nns)

   stat = file.getStat()
   cross=stat.getCrossSectionAccumulated()
   erro=stat.getCrossSectionErrorAccumulated();
   file.close()

print "Total Nr =",h1.allEntries()
c1 = HPlot("=HepSim=",600,700, 1,3) # plot histogram
c1.visible(True)
c1.setAutoRange()
#c1.setLogScale(1,True)
#c1.setRangeX(0,6000)
c1.setMarginLeft(90)
c1.setNameX("dR")
c1.setNameY("Entries / total")
h1.scale(1.0/h1.allEntries())
c1.draw(h1)
# h1.toTable() 
c1.drawStatBox(h1)

c1.cd(1,2)
c1.setAutoRange()
c1.setRangeY(0,0.3)
c1.setRangeX(0,30)
c1.setMarginLeft(90)
c1.setNameX("N(charged)")
c1.setNameY("Entries / total")
h2.scale(1.0/h2.allEntries())
c1.draw(h2)
c1.drawStatBox(h2)

c1.cd(1,3)
c1.setAutoRange()
c1.setRangeX(0,30)
c1.setMarginLeft(90)
c1.setNameX("N(stable)")
c1.setNameY("Entries / total")
h3.scale(1.0/h3.allEntries())
c1.draw(h3)
c1.drawStatBox(h3)

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