dijetmass_chah_hb.py (raw text file)# Validation script of HepMC: https://atlaswww.hep.anl.gov/hepsim/
# Run inside Jas4pp
# S.Chekanov (ANL)
from java.awt import Color,Font
from java.lang import *
from proto import FileMC
from jhplot import HPlot,P1D,HLabel,H1D
from jhplot.utils import FileList
from hepsim import HepSim
from java.util import ArrayList
from hephysics.hepsim import PromcUtil
from hephysics.jet import JetN2
from hephysics.particle import LParticle
import math,sys
ktjet=JetN2(0.4,"antikt",20)
dataset="tev13pp_mg5_chaH4FNS_sys"
tag="" # only truth level
url=""
NMax=0
if len(sys.argv)>1:
flist=FileList.get(sys.argv[1],"promc")
if (len(sys.argv)>2): NMax=int(sys.argv[2])
else:
sites=HepSim.urlRedirector(dataset)
url=sites[0]+"/"+tag
flist=HepSim.getList(url)
if len(flist)==0: print "Error: No input file!"; sys.exit(0)
else: print "Reading "+str(len(flist))+" files. Nr of files= ",NMax
print flist
# create file list using mass ranges
filelist=[]
histo=[] # m_jj mass
crosssection=[] # keeps cross sections
colors=[Color.black,Color.red, Color.blue, Color.green,Color.yellow]
mjjBins = [99,112,125,138,151,164,177,190, 203, 216, 229, 243, 257, 272, 287, 303, 319, 335, 352, 369, 387, 405, 424, 443, 462, 482, 502, 523, 544, 566, 588, 611, 634, 657, 681, 705, 730, 755, 781, 807, 834, 861, 889, 917, 946, 976, 1006, 1037, 1068, 1100, 1133, 1166, 1200, 1234, 1269, 1305, 1341, 1378, 1416, 1454, 1493, 1533, 1573, 1614, 1656, 1698, 1741, 1785, 1830, 1875, 1921, 1968, 2016, 2065, 2114, 2164, 2215, 2267, 2320, 2374, 2429, 2485, 2542, 2600, 2659, 2719, 2780, 2842, 2905, 2969, 3034, 3100, 3167, 3235, 3305, 3376, 3448, 3521, 3596, 3672, 3749, 3827, 3907, 3988, 4070, 4154, 4239, 4326, 4414, 4504, 4595, 4688, 4782, 4878, 4975, 5074, 5175, 5277, 5381, 5487, 5595, 5705, 5817, 5931, 6047, 6165, 6285, 6407, 6531, 6658, 6787, 6918, 7052, 7188, 7326, 7467, 7610, 7756, 7904, 8055, 8208, 8364, 8523, 8685, 8850, 9019, 9191, 9366, 9544, 9726, 9911, 10100, 10292, 10488, 10688, 10892, 11100, 11312, 11528, 11748, 11972, 12200, 12432, 12669, 12910, 13156]
nev=0
n=0
for j in flist:
if j.endswith("promc"):
filelist.append(j)
hh=H1D("sample="+str(n),mjjBins)
hh.setColor(colors[n])
histo.append( hh )
n=n+1
labels={}
print "Total number of mass ranges=",len(filelist)
Total=len(filelist)
for ipt in range(Total):
ptlist=filelist[ipt]
print ptlist
if (ptlist.find("tb_LO_CTEQ6L.promc")>-1): labels[ipt]="LO CTEQ6"
if (ptlist.find("tb_LO_NNPDF23.promc")>-1): labels[ipt]="LO NNPDF23"
if (ptlist.find("tb_NLO_NNPDF23.promc")>-1): labels[ipt]="NLO NNPDF23"
if (ptlist.find("tb_LO_NNPDF23_madspin.promc")>-1): labels[ipt]="LO NNPDF23 (MadSpin)"
file=FileMC(url+ptlist) # get input file
header = file.getHeader()
un=float(header.getMomentumUnit()) # conversion units
lunit=float(header.getLengthUnit())
print "Read=",url+ptlist
# if (nev>TotalEvents): print "Max Nr of events are done"; break # stop after some events
for i in range(file.size()):
nev=nev+1
if (nev%200==0):
print "Event=",nev
eve = file.read(i)
pa = eve.getParticles() # particle information
ve = eve.getEvent() # event information
leptons=[]
allpart=[]
# isolaton
for j in range(pa.getPxCount()):
if pa.getStatus(j)==1:
if (abs(pa.getPdgId(j))==12 or abs(pa.getPdgId(j))==14 or abs(pa.getPdgId(j))==15): continue # neutrino
par=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un,0)
allpart.append(par)
if (abs(pa.getPdgId(j))==11 or abs(pa.getPdgId(j))==13 ):
if (par.perp()>60 and abs(par.pseudoRapidity())<2.47): leptons.append(par)
# create isolated leptons
leptons_iso=[]
for lep in leptons:
phiL=lep.phi()
etaL=lep.rapidity()
ppT=0
for pp in allpart:
phi=pp.phi()
eta=pp.rapidity()
dR=math.sqrt( (phi-phiL)*(phi-phiL) + (eta-etaL)*(eta-etaL))
if (dR<0.2): ppT=ppT+pp.perp()
if (ppT/lep.perp() < 1.1): leptons_iso.append(lep) # 90% isolation
if len(leptons_iso)<1: continue
particles=PromcUtil.getParticleDList(file.getHeader(), pa, 1, 0, 1000);
ktjet.buildJets(particles);
jets=ktjet.getJetsSorted()
finaljets=[]
for jet in jets:
phi=jet.getPhi()
eta=jet.getRapidity()
isRemove=False
for lep in leptons_iso:
phiL=lep.phi()
etaL=lep.rapidity()
dR=math.sqrt( (phi-phiL)*(phi-phiL) + (eta-etaL)*(eta-etaL))
if (dR<0.2): isRemove=True # double counted
if (isRemove==False):
finaljets.append(jet)
# ktjet.printJets();
if (len(finaljets)>1): # mass of 2 leading jets
ptlead1=finaljets[0].perp()
ptlead2=finaljets[1].perp()
#print "jet1=",ptlead1
#print "jet2=",ptlead2
finaljets[0].add(finaljets[1])
mass=finaljets[0].mass()
histo[ipt].fill(mass)
# stat = file.getStat()
# xcross=stat.getCrossSectionAccumulated()
# cross=cross+xcross;
# erro=stat.getCrossSectionErrorAccumulated();
# file.close()
# cross=cross/Nfiles
# lumi=nev/cross
# crosssection.append( histo[ipt].entries()/lumi )
# histo[ipt].scale(1.0/lumi) # for cross section
# histo[ipt]=histo[ipt].getDividedByBinWidth() # get differential cross section
# print ipt, ") Cross section = %.3e pb"%cross, " Lumi=%.3e pb"%lumi
c1 = HPlot("HepSim",600,400) # plot histogram
c1.setNameX("M_{jj} [GeV]")
# c1.setNameY("d σ / d M_{jj} [pb / GeV]")
c1.setNameY("Entries / GeV")
c1.visible(True)
c1.setAutoRange()
c1.setRangeX(100,2200)
c1.setRangeY(0,500)
c1.setMarginLeft(100)
# c1.setLogScale(1,True)
# plot for mass=4 TeV
print labels
print "Cross sections for masses [pb]"
for ipt in range(Total):
print ipt,labels[ipt]
histo[ipt].setTitle(labels[ipt])
c1.draw(histo)
# h1.toTable()
#l1=HLabel("Pythia8 QCD σ=%.3e pb"%crossD, 0.5, 0.65, "NDC")
#l1.setFont(Font("Helvetica", Font.PLAIN, 14))
#c1.add(l1)
l2=HLabel("=HepSim=",0.75,0.84, "NDC")
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
l2.setColor(Color.gray)
c1.add(l2)
# create file/image using name of the file
name="output"
if len(sys.argv[0])>0: name=sys.argv[0].replace(".py","")
from jhplot.io import HBook
file=HBook(name+".jdat","w"); print name+".jdat created"
for i in histo: file.write(i)
file.close()
c1.export(name+".svg"); print name+".svg created"
# h1.toFile(name+".txt"); print name+".txt created"
# sys.exit(0)