higgs_0higgshiggs.py (raw text file)

# HepSim script using DataMelt http://jwork.org/dmelt
# Part of =HepMC= : https://atlaswww.hep.anl.gov/hepsim/
# S.Chekanov (ANL)
from java.lang import *
from proto import FileMC   # import FileMC
from jhplot.utils import FileList
from java.awt import Color,Font
from jhplot import  HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from hephysics.particle import LParticle 
from hepsim import HepSim
import sys, math 

"""
eta-> 2 Higgs -> 4 taus.
Calculate difference in pseudorapidity for H+H -> 4 taus
Provided by A.Kotwal
"""


url="";  TotalEvents=10000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/s0higgshiggs_alltau/"
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("Δ η (H0 to 4 τ)",50,-7,7)  
h1.setFill(True)
#h1.doc()               # check its methods

decays=[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
   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
      n=0
      higgs=[]    # keeps Higgs Lorentz vector 
      taus=[]
      for j in range(pa.getPxCount()):
           if (pa.getPdgId(j)==25):
                  higgs.append(LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un))
      if (len(higgs)>0):
             id1=len(higgs)-1;
             id2=len(higgs)-2;
             h1.fill(higgs[id1].pseudoRapidity()-higgs[id2].pseudoRapidity()) 

      # look at taus
      for j in range(pa.getPxCount()):
           if abs(pa.getPdgId(j))==15:
             dd=pa.getMother1(j)
             if (dd3:
             decays[0]= decays[0] +1 
   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,True)
c1.setMarginLeft(90)
c1.setNameX("Δ η (H0 to 4 τ)")
c1.setNameY("d σ / d Δ η [pb]")
c1.visible(True)
c1.setRange(-7,7,0.000000011,0.00001)

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 Higgs σ=%.3e pb"%cross, 0.3, 0.75, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 14))
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"

print "Branching 2H -> 4 tau      :", float(decays[0])/TotalEvents, " +- ", math.sqrt(float(decays[0]))/TotalEvents 
 

# xsec.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)