pythia_darkpi_pt.py (raw text file)

# HepSim script for hs-tools or Jas4pp 
# Look at pT of dark pions with M=1000 GeV
# Part of =HepMC= : https://atlaswww.hep.anl.gov/asc/hepsim/
# S.Chekanov (ANL)
from java.awt import Color,Font
from java.lang import *
from proto import FileMC            
from jhplot import HPlot,H1D,H2D,P1D,HLabel  # graphics
from jhplot.utils import FileList
from hephysics.particle import LParticle 
from hephysics.hepsim import PromcUtil
from hephysics.jet import JetN2 
from hepsim import HepSim
import sys,math 

TotalEvents=10000


dataset="tev27pp_pythia8_darkpion_scan"
tag="" # only truth level
url=""
NMax=0
if len(sys.argv)>1:
   flist=FileList.get(sys.argv[1],"promc")
   if (len(sys.argv)>2): NMax=int(sys.argv[2])
else:
   sites=HepSim.urlRedirector(dataset)
   url=sites[0]+"/"+tag
   flist=HepSim.getList(url)
if len(flist)==0: print "Error: No input file!"; sys.exit(0)
else: print "Reading "+str(len(flist))+" files. Nr of files= ",NMax


#h1.doc()               # check its methods
h1= H1D("pT [GeV]",100,0,2000)


c1 = HPlot("=HepSim=",500,500)
c1.visible(True)
c1.setAutoRange(True)
c1.setMarginLeft(90)
c1.setLegend(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
   print "Read=",flist[m]
   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%1000==0):
           if (Nfiles==1): print "Event=",nev
           else: 
                print "Event=",nev," done=",int((100.0*m)/Nfiles),"%"
                c1.clearData()
                #c1.setNameX("p_{T}(lead. jet) [GeV]")
                #c1.setNameY("Entries")
                c1.draw(h1)
                print "Update canvas"
      eve = file.read(i)
      pa = eve.getParticles()    # particle information
      nt=0; xsum=0; 
      darkpi={}
      #print " "
      for j in range(pa.getPxCount()):
        apdg=abs(pa.getPdgId(j))
        if apdg<4:   # look at quarks. Can use  mesons too to  estimate ctau 
           ida=pa.getMother1(j)
           if (abs(pa.getPdgId(ida))==4900111 or abs(pa.getPdgId(ida))==4900113):
             #print apdg # pa.getPdgId(ida)
             xL=pa.getX(j)/lunit
             yL=pa.getY(j)/lunit
             zL=pa.getZ(j)/lunit
             ctau2=xL*xL+yL*yL+zL*zL 
             ctau=math.sqrt(ctau2)
             par=LParticle(pa.getPx(ida)/un,pa.getPy(ida)/un,pa.getPz(ida)/un,pa.getEnergy(ida)/un,pa.getMass(ida)/un)
             nt=nt+1
             h1.fill(par.perp()) # calculate transverse mementa 
   stat = file.getStat()
   cross=stat.getCrossSectionAccumulated()
   erro=stat.getCrossSectionErrorAccumulated();
   file.close()

# show final cross section
c1.clearData()
c1.clear()
#c1.setLogScale(1,1)
#c1.setRangeY(0,01)
c1.setRangeX(0,1000)
c1.setMarginLeft(90)
c1.setNameX("p_{T} [GeV]")
c1.setNameY("Events")

c1.draw(h1)
print h1.getStat()
#l1=HLabel("σ=%.3e pb"%cross, 0.4, 0.75, "NDC")
#l1=HLabel("=%.3e GeV"%h2.mean(), 0.4, 0.75, "NDC")
#l1.setFont(Font("Helvetica", Font.PLAIN, 15))
#c1.add(l1)

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

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