HepSim

Repository with Monte Carlo simulations for particle physics

higgsjet_gamgam_mcfm.py (raw text file)
# HepSim script using DataMelt http://jwork.org/dmelt
# Shows Higgs x-section with PDF uncertainties
# Higgs is built from 2 gammas
# Part of =HepMC= : https://atlaswww.hep.anl.gov/asc/hepsim/
# S.Chekanov (ANL)

from java.awt import Color,Font 
from jhplot import  HPlot,P1D,HLabel,H1D
from java.lang import *
from proto import FileNLO 
from jhplot.utils import FileList
from hephysics.particle import LParticle # physics
from hepsim import HepSim
import os,sys

url="";  TotalEvents=10000000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/higgsjet_gamgam_mcfm/"
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


bins=[0,30,60,90,120,150,180,210,240,280,320,360,400,500,600,700,800,1000,1500,2000]
minJetPt=50
nPDF=52     # 52 histograms with different PDF
histo=[]
for i in range(nPDF):
      h1=H1D("PT(H0) pdf="+str(i),bins)
      histo.append(h1)

#histo[0].doc()   # check its methods
nev=0;  Nfiles=len(flist)
for m in range(Nfiles):              # loop over all files in this list    
   file=FileNLO(url+flist[m])        # get input file
#  if (m>50): break                   # uncomment for debug
   header = file.getHeader()
   un=float(header.getMomentumUnit()) # convertion units
   lunit=float(header.getLengthUnit())
   nIterations=header.getWeight()     # MCFM Nr of iterations
   NXF=nIterations*Nfiles
   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%20000==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
      gweight=ve.getFdata(0)     # global weight. 1-gg, 2-gq, 3-qq, 4-qqb
      weight=gweight/NXF         # take into account iterations & Nr of files
      gammas=[];  jets=[] 
      for j in range(pa.getPxCount()):
           p=LParticle(pa.getPx(j)/un, pa.getPy(j)/un, pa.getPz(j)/un, pa.getEnergy(j)/un)
           if (pa.getPdgId(j)==22): gammas.append(p) 
           if (pa.getPdgId(j)==21 and p.perp()>minJetPt): jets.append(p) 
      mass=0
      higgs=gammas[0]
      if (len(gammas)>1):
                   higgs.add(gammas[1])
                   mass=higgs.calcMass()
      if (mass>120 and mass0):  # ask for a jet above 50 GeV  
                     for k in range(nPDF):
                           delta = 1.0-(ve.getIdata(k)/1000.) # PDF errors (1-PDF(i)/PDF(0))*1000
                           histo[k].fill(higgs.perp(),weight*delta)

   file.close()
 
cross=histo[0].sumAllBinHeights()/1000 # convert to pb
print "Total cross section (pb)=",cross
data=[]
for k in range(nPDF):
    newD=histo[k].getDividedByBinWidth()
    newD.scale(0.001) # convert to pb 
    data.append(P1D(newD))

# create P1D with errors evaluated from an array of P1D
xsec=P1D("MCFM CT10")
for j in range(data[0].size()):
     x=data[0].getX(j)
     y=data[0].getY(j)
     errX=0.5*(bins[j+1]-bins[j])
     up=0; dw=0;
     for n in range(1,nPDF):
        d=data[n].getY(j)-y
        if (d>0): up=up+d*d
        if (d0: name=sys.argv[0].replace(".py","")
file=HBook(name+".jdat","w"); print name+".jdat created"
file.write(xsec)
file.write(errors)
file.close()
c1.export(name+".svg");    print name+".svg created"
# xsec.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)

HEP.ANL.GOV