zboson_ee_mcfm.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.awt import Color,Font
from jhplot import  HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from java.lang import *
from proto import FileNLO # import FileNLO
from jhplot.utils import FileList
from hephysics.particle import LParticle
from hepsim import HepSim
import sys

url="";  TotalEvents=100000000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/zboson_ee_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=[10,20,30,40,50,60,70,80,90,100,120,140,160,200,240,280,320,360,400,500,600,700,800,1000,1500,2000]
nPDF=52
histo=[]
for i in range(nPDF):
      h1=H1D("pdf="+str(i),bins)
      histo.append(h1)

nev=0;  Nfiles=len(flist)
for m in range(Nfiles):               # loop over all files in this list    
#   if (m>50): break                  # debugging
   file=FileNLO(url+flist[m])         # get input file
   header = file.getHeader()
   un=float(header.getMomentumUnit()) # GeV 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
      weight=ve.getFdata(0)      # global weight. 1-gg, 2-gq, 3-qq, 4-qqb
      for j in range(ve.getIdataCount()): 
                pdferrors= ve.getIdata(j)      # PDF errors (1-PDF(i)/PDF(0))*1000 
      weight=weight/NXF            # take into account iterations & Nr of files
      ele=[]
      for j in range(pa.getPxCount()):
         if (Math.abs(pa.getPdgId(j))==11):
           p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un,0)
           ele.append(p)
      if (len(ele)>1):
           zboson=ele[0]
           zboson.add(ele[1])
           mass=zboson.calcMass()
           if (mass>80 and mass<100): 
             for k in range(nPDF):
                     delta= 1.0-(ve.getIdata(k)/1000.) # PDF errors (1-PDF(i)/PDF(0))*1000
                     histo[k].fill(zboson.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 (d<0): dw=dw+d*d
     xsec.add(x,y,errX,errX,Math.sqrt(up),Math.sqrt(dw))
xsec.setErr(1)       # show errors
xsec.setColor(Color.blue)


errors=P1D("Relative errors")
for j in range(xsec.size()):
   y=xsec.getY(j)
   if y !=0:
     errors.add(xsec.getX(j),1,xsec.getXleft(j),xsec.getXright(j),xsec.getYupper(j)/y,xsec.getYlower(j)/y)
errors.setErr(1)
# errors.setErrFillColor(Color.blue,0.3)

c1 =HPlot("=HepSim=",400,600, 1,2)
c1.visible(True)
c1.resizePad(0,0,1.0, 1.5)
c1.resizePad(0,1,1.0, 0.5)

xmin=0
xmax=2000
c1.cd(1,1)
c1.setNameY("d σ / d p_{T}(Z) [pb/GeV]")
c1.setNameYpos(-0.01, 0.93)
c1.setMarginLeft(55)
c1.setMarginBottom(0.0)
c1.setTicLabels(0, False)
c1.setRange(xmin,xmax,0.0000011,1000)
c1.setLogScale(1,1)

tot=0
for i in range(xsec.size()):
  tot=tot+xsec.getY(i)*(bins[i+1]-bins[i])
l1=HLabel("σ(Z(ee))=%.3e pb"%cross, 0.15, 0.7, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 15))
c1.add(l1)

l2=HLabel("=HepSim=",0.5,0.8, "NDC")
l2.setColor(Color.gray)
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l2)
c1.draw(xsec)


c1.cd(1,2)
c1.setRange(xmin,xmax,0.7,1.3)
c1.setLegend(False)
c1.setNameY("Ratio")
c1.setNameYpos(-0.01, 0.9)
c1.setMarginTop(0.0)
c1.setMarginLeft(55)
c1.setNameX("p_{T}(Z) [GeV]")
c1.setNameXpos(0.55, 0.1)


from jhplot.shapes  import Line
line = Line(xmin,1.0, xmax, 1.0)
line.setPosCoord("USER")
line.setTransparency(0.5)
c1.add(line)
c1.draw(errors)

# 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.write(errors)
file.close()
c1.export(name+".svg");    print name+".svg created"
# xsec.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)

• Get involved
• Full Search
• Experiments
  • LHC 8 TeV (pp)
  • LHC 13 TeV (pp)
  • HL-LHC (pp)
  • HE-LHC (pp)
  • FCC-ee (e+e-)
  • FCC-hh (pp)
  • SppC (pp)
  • CEPC (e+e-)
  • CLIC (e+e-)
  • TESLA (e+e-)
  • ILC (e+e-)
  • Muon collider
  • HERA (ep)
  • EIC (ep)
• Manual
• Mirrors
  • NERSC.GOV
  • CERN.CH
  • JLAB.ORG
• Tools
  • Detectors
  • RecoTags
  • EVGN validation
  • RECO validation
  • Jas4PP
• About
  • About
  • File statistics
• Login

HEP.ANL.GOV