higgs_h2all_pythia8.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 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 

"""
Double Higgs production in a BSM scenario.
Calculate difference in pseudorapidity and estimate branching ratios to 4 b, 2b+2tau, 4 tau
"""


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

decays=[0,0,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 
      higgs_id=[] # keeps ID of decays 

      for j in range(pa.getPxCount()):
           if (pa.getPdgId(j)==25):
              id=pa.getDaughter1(j)
              id_m=pa.getPdgId(id)
              if (id_m != 25):  
                      higgs.append(LParticle(pa.getPx(id)/un,pa.getPy(id)/un,pa.getPz(id)/un,pa.getEnergy(id)/un))
                      higgs_id.append([abs(pa.getPdgId(pa.getDaughter1(j))),  abs(pa.getPdgId(pa.getDaughter2(j)))])  
      if (len(higgs)==2):
             h1.fill(higgs[0].pseudoRapidity()-higgs[1].pseudoRapidity()) 
             hig1=higgs_id[0]
             hig2=higgs_id[1]
             # print hig1[0],hig1[1],hig2[0],hig2[1] # print decays ID 
             if (hig1[0]==5 and hig1[1]==5) and (hig2[0]==5 and hig2[1]==5): decays[0]=decays[0]+1 # 4 b decays 
             if ( (hig1[0]==5 and hig1[1]==5 and hig2[0]==15 and hig2[1]==15) or  (hig1[0]==15 and hig1[1]==15 and hig2[0]==5 and hig2[1]==5) ):  decays[1]=decays[1]+1 # 2 b + 2 tau decays 
             if (hig1[0]==15 and hig1[1]==15 and hig2[0]==15 and hig2[1]==15): decays[2]=decays[2]+1 # 4 tau decays 

    
             # if (pa.getPdgId(m_id) == 25): print "Higgs" 
             #p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
             #h1.fill(p.perp()) 
             #n=n+1
      # print "Nr of higgs=",n 
   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)")
c1.setNameY("d σ / d Δ η [pb]")
c1.visible(True)
c1.setRange(-5,5,0.00011,10)

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 -> 2 bbar     :", float(decays[0])/TotalEvents, " +- ", math.sqrt(float(decays[0]))/TotalEvents 
print "Branching 2H -> 2 bbar+2tau:", float(decays[1])/TotalEvents, " +- ", math.sqrt(float(decays[1]))/TotalEvents
print "Branching 2H -> 4 tau      :", float(decays[2])/TotalEvents, " +- ", math.sqrt(float(decays[2]))/TotalEvents 
 

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