HepSim

Repository with Monte Carlo simulations for particle physics

gamma_jetphox_ptbins.py (raw text file)
# HepSim script using DataMelt http://jwork.org/dmelt
# Shows gamma cross section with 52 PDF uncertainties
# 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  
from jhplot.utils import FileList
from hephysics.particle import LParticle
import os, sys
import glob


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); 
           


flist=[] # file list
ptbins=[[200,1500]]
url=None
print "Nr of arguments=",len(sys.argv)
if len(sys.argv)>1: # if one argument is given, read all files from a directory 
   DIR=sys.argv[1]
   print DIR
   flist200=glob.glob(DIR+"*_mu1_200_600_*promc")
   flist600=glob.glob(DIR+"*_mu1_600_1200_*promc")
   flist1200=glob.glob(DIR+"*_mu1_1200_2400_*promc")
   flist2400=glob.glob(DIR+"*_mu1_2400_4800_*promc")
   flist4800=glob.glob(DIR+"*_mu1_4800_15000_*.promc")
   print len(flist200), len(flist600), len(flist1200), len(flist2400), len( flist4800)
   if (len(flist200)*len(flist600)*len(flist1200)*len(flist2400)*len(flist4800)==0):
               print "Error: some files are missing !"; sys.exit(0)
   flist= [flist200, flist600, flist1200, flist2400, flist4800]
   ptbins=[[200,600],[600,1200],[1200,2400],[2400,4800],[4800,15000]]
else:
   url="http://mc.hep.anl.gov/asc/hepsim/events"
   flist=[url+"/pp/100tev/gamma_jetphox/ggd_mu1_1000_6000_run0_atlas50.promc"]
   print "Reading 1 file from WWW=",flist[0]

isoCone=0.4
isoET=7
EtaMax=2.5

nPDF=41
bins=[200,300,400,500,600,700,800,900,1000,1200,1400,1600,1800,2000,2200,2400,2600,3000,3400,3800,4200,4800,5000,6000,7000,8000,10000,14000]
pthisto=[]

for ipt in range(len(flist)):
 ptlist=flist[ipt]
 histoD=[] # histograms with different PDF for direc 
 histoF=[] # histograms with different PDF for fragm 
 print "Dealing with="+str(ptbins[ipt])+" GeV" 
 for i in range(nPDF):
      h1=H1D("PT(direct) pdf="+str(i),bins)
      histoD.append(h1)
      h2=H1D("PT(fragm) pdf="+str(i),bins)
      histoF.append(h2)

 nfilesF=0; nfilesD=0; 
 nev=0;  Nfiles=len(ptlist)
 for m in range(len(ptlist)):              # loop over all files in this list    
#   if (m>10): break;
   bsize=os.path.getsize(ptlist[m])
   if (bsize-1:
                    fragm=True
                    nfilesF=nfilesF+1
   else: nfilesD=nfilesD+1 
   #print "nfilesD=",nfilesD, " nfilesF=",nfilesF
   header = file.getHeader()
   un=float(header.getMomentumUnit()) # GeV units
   lunit=float(header.getLengthUnit())
   stat = file.getStat()
   cross=stat.getCrossSectionAccumulated()
   nevents=stat.getNTried()
   lumi=nevents/cross
   norm=1.0/lumi
   if m==0 and ipt==0:
      print "ProMC v=",file.getVersion(), "M unit=",un,"L unit=",lunit
   for i in range(file.size()):
      nev=nev+1  
      if (nev%20000==0): 
           print "Event=",nev," done=",int((100.0*m)/Nfiles),"%"
      eve = file.read(i) 
      pa = eve.getParticles()    # particle information
      ve = eve.getEvent()        # event information
      wi=ve.getFdata(0)
      if (ve.getIdataCount() != nPDF): print "Error in PDF set!" 

      for j in range(pa.getPxCount()):
         if (pa.getPdgId(j)==22):
           p=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un)
           pt=p.perp()
           eta=p.pseudoRapidity()
           if (pt>ptbins[ipt][0] and pt0) :
   crossF=histoF[0].sumAllBinHeights()/nfilesF

 print "Direct cont =",crossD," from ",nfilesD," files" 
 print "Fragm  cont =",crossF," from ",nfilesF," files"    

 data=[]
 for k in range(nPDF):
    if (nfilesD>0) :
      newD=histoD[k].getDividedByBinWidth()
      newD.scale(1.0/nfilesD)
      newD.scaleErrors(0)
    if (nfilesF>0) :
      newF=histoF[k].getDividedByBinWidth()
      newF.scale(1.0/nfilesF)
      newF.scaleErrors(0)
    if nfilesD>0 and nfilesF>0:
       hh=newD.oper(newF,"+")
       data.append(P1D(hh))
    if nfilesD>0 and nfilesF==0:
       data.append(P1D(newD))

 if (len(data)==0):
        print "Error: No files!"; sys.exit(0)
 
# create P1D with 2nd level errors evaluated from an array of P1D
 #xsec=data[0]
 #xsec=xsec.getSys(data)
 #xsec.setErrToZero(1) # set all stat. error to 0
 xsec=P1D("Cross sections:"+str(ptbins[ipt]))
 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       # scale+PDF 
     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","")
## save in a few formats 
file=HBook(name+".jdat","w"); print name+".jdat created"
file.write(pp)
file.write(errors)
for k in range(len(flist)): file.write(pthisto[k])
file.close()
c1.export(name+".svg");    print name+".svg created"
# xsec.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)

HEP.ANL.GOV