higgs_gammaphi.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 
from jhplot.utils import FileList
from hephysics.particle import LParticle 
from hepsim import HepSim
import  sys

url="";  TotalEvents=10000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/14tev/higgs/pythia_gammaphi/"
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("Mass(γ-φ)",200,110,140)   # create a histogram
h1.setFill(True)
#h1.doc()               # check its methods

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
      gammas=[]
      kP=[]
      kM=[]
      for j in range(pa.getPxCount()): 
        if (pa.getStatus(j)==1):
           if (pa.getPdgId(j)==22): 
                p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
                if (p.perp()>10): gammas.append(p)
           if (pa.getPdgId(j)==-321):
                p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
                if (p.perp()>1): kM.append(p)
           if (pa.getPdgId(j)==321):
                p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
                if (p.perp()>1): kP.append(p)

      phimesons=[] # make phi from K+ and K- in some mass window 
      for k1 in kP:
               for k2 in kM:
                  k1.add(k2)
                  mass=k1.calcMass()
                  if (mass>0.95 and mass<1.1): phimesons.append(k1)
      # print len(phimesons)

      for p in phimesons:
               for g in gammas:
                 p.add(g)
                 mass=p.calcMass()
                 h1.fill(mass)
                 # print mass  

   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
print "Total Nr of Higgs=",h1.allEntries()
h1.scale(1.0/(h1.getBinSize()*lumi))
c1 = HPlot("=HepSim=",600,400) # plot histogram
c1.setNameX("M [GeV]")
c1.setNameY("d σ / d M [pb / GeV]")
c1.visible(True)
c1.setAutoRange()
#c1.setRangeX(0,200)
#c1.setRangeY(0,1000000)
c1.setMarginLeft(70)

c1.draw(h1)
# h1.toTable() 

#l1=HLabel("Pythia8 σ=%.3e pb"%cross, 0.4, 0.8, "NDC")
#c1.add(l1)
l2=HLabel("=HepSim=",0.75,0.86, "NDC")
l2.setFont(Font("Helvetica", Font.PLAIN, 12))
l2.setColor(Color.gray)
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(h1)
file.close()
c1.export(name+".svg");    print name+".svg created"
# h1.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)