HepSim

Repository with Monte Carlo simulations for particle physics

gamma_mcfm_ptscale_alpha020.py (raw text file)
# Shows MCFM cross section for inclusive gamma with 52 PDF uncertainties
# Part of HepSim:  http://atlaswww.hep.anl.gov/hepsim/
# Runs in DataMelt http://jwork.org/dmelt/
# S.Chekanov (ANL)
from java.lang import *
from proto import FileNLO # import FileNLO
from java.awt import Color,Font
from jhplot import  HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from jhplot.utils import FileList
from hephysics.particle import LParticle
from hepsim import HepSim
import math 
import sys

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="";  TotalFiles=-1
# set TotalFiles to negative to run over all
www="http://mc.hep.anl.gov/asc/hepsim/events/pp/100tev/gamma_mcfm_ptscale/"
samples=["20","35","55","75","105","150","200","300","400","550","750","1100","1500","2000","2500","3000","4000","5000","6000","7000","8000","9000","10000","15000"]


data=[]

import glob

for i in range(len(samples)-1):
    data.append(glob.glob('./data/alpha020/pt'+samples[i]+'_*.promc') )

if len(data)==0: print "Error: No input file!"; sys.exit(0)
else: print "Reading "+str(len(data))+" samples. Nr = ",len(data)


# isolation cut as in atlas
def getIsoCut(pt):
    return 4.8+4.2E-03*pt;

# define Pi
Pi=math.pi
Pi2=2*math.pi

# 52 histograms with different PDF (CT10)  
nPDF=1   
#histoEta1= # Eta 0):  Nfiles=2

 for m in range(Nfiles):             # loop over all files in this list    
   file=FileNLO(url+flist[m])        # get input file
   header = file.getHeader()
   #   print header
   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
   for i in range(file.size()):
      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)/NXF  # global weight. 1-gg, 2-gq, 3-qq, 4-qqb 
      pdf=weight
      #for k in range(nPDF):
      #             delta=1.0-(ve.getIdata(k)/1000.) # PDF errors (1-PDF(i)/PDF(0))*1000
      #             pdf.append(weight*delta)

      gammas=[]
      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)
           gammas.append(p) 
 
      for p in gammas: 
       esum=0
       pt_gam=p.perp()
       eta_gam=p.pseudoRapidity()
       abseta_gam=Math.abs(eta_gam)
       phi_gam=p.phi()
       if pt_gam>float(samples[ii]) and pt_gamPi): adphi=Pi2-adphi
           ddr=Math.sqrt(deta*deta+adphi*adphi);
           # print ddr 
           if (ddr0: print "Isolation=",esum 
        if (esum

HEP.ANL.GOV