truth_q2_1gev_lcio.py (raw text file)

# Validation script: Read LCIO files
# Use Jas4pp as: fpad script.py
# Task: Extract truth-level Q2 variable and compare with Q2 from PFA 
# S.Chekanov (ANL)

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 java.awt import Color,Font
from jhplot import *  # import graphics
from jhplot.io import HBook
from hephysics.particle import LParticle
from java.lang import System
from java.io import *
from java.net import URL
from hepsim import HepSim
import sys,math 
from  java.awt import Color
from org.apache.commons.io import FileUtils

NFilesMax=10 # max number of files 
dataset="gev35ep_lepto6_dis1q2"
tag="rfull055"

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("Q2(ele true)",100,0,100)
h2=H1D("Q2(ele reco)",100,0,100)

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
    if (nEvent%100==0): print "# Event: ",nEvent
    params=evt.getParameters()
    Q2true=params.getFloatVal("EVGEN:DIS:Q2")  
    h1.fill(Q2true);

    # get electron beam energy
    col = evt.getCollection("MCParticle")
    nMc=col.getNumberOfElements()
    e_beam=0;
    for i in range(nMc): # loop over all particles 
      par=col.getElementAt(i)
      momentum = par.getMomentum()
      ee=par.getEnergy()
      if (abs(par.getPDG())==11):  
                            e_beam=ee;
                            break
    # get reconstructed electrons using PFA
    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()
          p=LParticle(p4[0],p4[1],p4[2],ee)
          p.setCharge(charge)
          if (abs(par.getPDG())==11): 
                                particles.append(p)  
                                break
    if (len(particles)>0):
          ele=particles[0]       
          mag=ele.mag()
          theta=math.atan(ele.perp()/p.pz())
          if(theta<0): theta+=3.14159
          Q2ele=2.*e_beam*ele.e()*(1+math.cos(theta))
          # print Q2true,Q2ele,e_beam 
          h2.fill(Q2ele);
 
    del col
    del evt 
  reader.close() # close the file
  del reader
  xfile.delete() 
  System.gc() # if you low on memory, manual JVM cleanup for every 5th file 

c1=HPlot() # plotting part
c1.visible()
c1.setAutoRange()
c1.setMarginLeft(90)
c1.setNameX("Q^{2} [GeV^{2}]")
c1.setNameY("Events")
c1.setRangeX(0,10)
c1.draw(h1)
p2=P1D(h2)
p2.setColor(Color.blue)
c1.draw(p2) # show reco Q2 

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)