HepSim

Repository with Monte Carlo simulations for particle physics

jets_nlojetpp.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 FileNLO  # import FileNLO
from jhplot.utils import FileList
from hephysics.particle import LParticle
import glob,os, sys
from hepsim import HepSim

from javax.swing import JOptionPane
if JOptionPane.showConfirmDialog(None, "This macro works for local files. Quit?","Warning", JOptionPane.OK_CANCEL_OPTION) == JOptionPane.OK_OPTION: System.exit(0);



url=""
if len(sys.argv)>1: # if one argument is given, read all files from a directory 
   DIR=sys.argv[1]
   flist50=glob.glob(DIR+"*_pt50_*promc")
   flist100=glob.glob(DIR+"*_pt100_*promc")
   flist200=glob.glob(DIR+"*_pt200_*promc")
   flist400=glob.glob(DIR+"*_pt400_*promc")
   flist800=glob.glob(DIR+"*_pt800_*promc")
   flist1600=glob.glob(DIR+"*_pt1600_*promc")
   flist3200=glob.glob(DIR+"*_pt3200_*promc")
   flist6400=glob.glob(DIR+"*_pt6400_*promc")
   if (len(flist50)*len(flist100)*len(flist200)*len(flist800)*len(flist1600)*len(flist3200)*len(flist6400)==0):
               print "Error: some files are missing !"; sys.exit(0)
   flist= [ flist50, flist100, flist200,  flist400, flist800,   flist1600,  flist3200, flist6400]
else:
   url="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/jets_nlojetpp/"
   flist=HepSim.getList(url)
if len(flist)==0: print "Error: No input file!"; sys.exit(0)

minPT=50
maxETA=3.5
print "Inc. cross section. minPT=",minPT," max Eta=",maxETA

SYST=1+7+40 # 0, central, 7- scale uncertainties, 40-PDF (MSTW) 
bins=[50,100,200,400,800,1600,3200,6400,12800]
# bins=[50,75,100,150,200,300,400,600,800,1200,1600,2400,3200,6400]
nev=0
ptbins=[[50,100], [100,200],[200,400], [400,800],[800,1600],[1600,3200],[3200,6400],[6400,10000]]
pthisto=[]
print "Process all pT ranges:",ptbins

for ipt in range(len(flist)):
  ptlist=flist[ipt]
  print "Dealing with="+str(ptbins[ipt])+" GeV"
  Nfiles=len(ptlist)
  totruns=0
  avcross=0
  nfiles=0
  histo=[]    # histograms with different PDF for direct
  for i in range(SYST):
      h1=H1D("PT pdf="+str(i),bins)
      histo.append(h1)

  for m in range(len(ptlist)):           # loop over all files in this list    
#    if (m>10): break #  debugging mode 
    bsize=os.path.getsize(ptlist[m])
    if (bsizeptbins[ipt][0] and pt

HEP.ANL.GOV