zpeak_cepc_lcio.py (raw text file)

# Plot Z peak from e+e- using PFA electrons 
# Use Jas4pp to run this script. 
# author: S.Chekanov (ANL)

from org.lcsim.event import *
from org.lcsim.util import *
from hep.physics.vec import *
from java.util import *
from java.io import *
from java.awt import Color,Font
from org.lcsim.lcio import LCIOReader
from hep.io.sio import SIOReader
from hep.lcio.implementation.sio import SIOLCReader
from hep.lcio.implementation.io import LCFactory
from hep.lcio.event import *
from hep.lcio.io import *
from jhplot import *  # import graphics
from hephysics.particle import LParticle 
from math import *
from java.net import URL
from hepsim import HepSim
from jhplot.io import HBook
import sys
from org.apache.commons.io import FileUtils

# make list of files..
NFilesMax=10 # max number of files 
dataset="gev250ee_pythia6_zpole_ee"
tag="rfull101"

www=""
if len(sys.argv)>1:
   flist=FileList.get(sys.argv[1],"slcio")
   if (len(sys.argv)>2): FilesMax=int(sys.argv[2])
else:
   sites=HepSim.urlRedirector(dataset)
   www=sites[0]+"/"+tag
   flist=HepSim.getList(www)
if len(flist)==0: print "Error: No input file!"; sys.exit(0)
else: print "Reading "+str(len(flist))+" files. Nr of files= ",NFilesMax


h1=H1D("Mass",70,50,120)
factory = LCFactory.getInstance()
nEvent=0
for n in range(len(flist)):
  if (n+1>NFilesMax): break
  print "Open file=", flist[n]
  url=URL(www+"/"+flist[n])
  xfile=File("/tmp/tmp.slcio");
  FileUtils.copyURLToFile(url, xfile)
  reader = factory.createLCReader()
  reader.open(xfile.getAbsolutePath())
  while(1):
     evt=reader.readNextEvent()
     if (evt == None): break
     nEvent=nEvent+1
     # print " file event: ",evt.getEventNumber(), " run=",evt.getRunNumber()
     if (nEvent%100==0): print "# Event: ",nEvent
     strVec = evt.getCollectionNames()
     if nEvent == 1:
            for col in  strVec:
                           print col           
     col = evt.getCollection("PandoraPFOCollection")
     nPFA = col.getNumberOfElements()

     particles=[]
     for i in range(nPFA):
          # http://www.lcsim.org/sites/lcsim/apidocs/org/lcsim/lcio/SIOReconstructedParticle.html
          # http://www.lcsim.org/sites/lcsim/apidocs/org/lcsim/event/ReconstructedParticle.html
          pa=col.getElementAt(i)
          charge=pa.getCharge()
          p4=pa.getMomentum()  
          ee=pa.getEnergy() 
          typep=pa.getType() # type
          p=LParticle(p4[0],p4[1],p4[2],ee)
          p.setCharge(charge)
          if (abs(typep)==11):  particles.append(p) # reconstructed electrons 
     if (len(particles)>1):
       for i1 in range(0,len(particles)-1):
          p1=particles[i1]
          charge1=p1.getCharge()
          for i2 in range(i1+1,len(particles)):
              p2=particles[i2]
              charge2=p2.getCharge()
              p1.add(p2)
              if (charge1*charge2<0):
                  # print p1.calcMass()
                  mass=p1.calcMass()
                  h1.fill(mass)

reader.close();


c1=HPlot()
c1.visible()
c1.setAutoRange()
c1.setMarginLeft(100)
c1.setNameX("Mass [GeV]")
c1.setNameY("Events")
c1.setRangeX(50,120)
c1.draw(h1)

l2=HLabel("=HepSim=",0.9,0.86, "NDC")
l2.setColor(Color.gray)
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l2)

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");    print name+".svg created"
# xsec.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)