pythia_Ggamma_15tev_mT.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 
import sys
from hepsim import HepSim

url="";  TotalFiles=-1
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/13tev/monophoton15tev/"
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): TotalFiles=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"

h1= H1D("Transverse mass",100,0,1000)  
h1.setFill(True)
#h1.doc()               # check its methods

cross=0; nev=0;  Nfiles=len(flist)
if (TotalFiles>-1): 
               Nfiles=TotalFiles
               print "Debug mode for ",TotalFiles," files"
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%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
      gamma=[]
      Graviton=[]
      px=0;
      py=0;
      for j in range(pa.getPxCount()): 
         if (pa.getStatus(j)==1):
           px=px+pa.getPx(j)/un
           px=py+pa.getPy(j)/un
           pt=Math.sqrt(px*px+py*py)
           if (pa.getPdgId(j)==22 and pt>10):   
              p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un,pa.getMass(j)/un)
              gamma.append(p) 
      if (len(gamma)>0): 
        pxG=-px-gamma[0].px() 
        pyG=-py-gamma[0].py() 
        ptG=Math.sqrt(pxG*pxG+pyG*pyG)
        phiG=Math.atan2(pyG,pxG)
        ptGamma=gamma[0].perp()
        phiGamma=gamma[0].phi()
        mT=Math.sqrt(2*ptGamma*ptG*(1-Math.cos(phiGamma-phiG))) 
        h1.fill(mT)


   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 = HPlot("=HepSim=",500,500)
c1.setLogScale(1,0)
c1.setMarginLeft(80)
c1.setNameX("M_{T} [GeV]")
c1.setNameY("d σ / d p_{T} [pb/GeV]")
c1.visible(True)
c1.setRange(0,1000,0.0,0.05)

hNew=h1.getDividedByBinWidth()
hNew.scale(1.0/lumi)
xsec=P1D(hNew)
xsec.setErr(1)       # show errors
xsec.setColor(Color.blue)
c1.draw(xsec)
# h1.toTable() 

l1=HLabel("Pythia8 G+γ σ=%.3e pb"%cross, 0.3, 0.75, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 15))
c1.add(l1)

l2=HLabel("=HepSim=",0.6,0.87, "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","")
file=HBook(name+".jdat","w"); print name+".jdat created"
file.write(xsec)
file.close()
c1.export(name+".svg");    print name+".svg created"
# xsec.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)