pythia8_gluino_llp.py (raw text file)

# HepSim script using DataMelt http://jwork.org/dmelt
# Part of =HepMC= : https://atlaswww.hep.anl.gov/hepsim/
# S.Chekanov (ANL)
from java.awt import Color,Font
from java.lang import *
from java.util import ArrayList 
from proto import FileMC            
from jhplot import HPlot,H1D,P1D,HLabel 
from jhplot.utils import FileList
from hephysics.hepsim import PromcUtil
from hephysics.jet import ParticleD,JetN2 
from hepsim import HepSim
import math,sys

"""
Look at SUSY scenario of q g -> squark gluino, when gluino is long-lived (tau0=2 m).
Look at jets from LLP. 
"""

ktjet=JetN2(0.5,"antikt",100) # anti-KT,E-mode, R=0.5,min pT=100 GeV,fast
print ktjet.info()


url="";  TotalEvents=6000
default_www="http://mc.hep.anl.gov/asc/hepsim/events/pp/13tev/pythia8_gluino_llp/"
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


h1= H1D("PT (good jet)",50,0,2000)
h2= H1D("PT (from LLP)",50,0,2000)
h2.setColor(Color.red)
h3= H1D("Delta phi",100,0,5)

#h1.doc()               # check its methods

cross=0; nev=0;  Nfiles=len(flist)
for m in range(Nfiles):              # loop over all files in this list    
   file=FileMC(url+flist[m])         # get input file
   header = file.getHeader()
   un=float(header.getMomentumUnit()) # conversion units
   lunit=float(header.getLengthUnit())
   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%500==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

      particles=ArrayList()      # list to keep stable short lived particles 
      ll_particles=ArrayList()   # list to keep stable long-lived particles 
      for j in range(pa.getPxCount()):
          apdg=abs(pa.getPdgId(j))
          if (apdg==12 or apdg==14 or apdg==16): continue # no neutrino
          if (apdg==13): continue                         # no muons 
          if (pa.getStatus(j)==1):
             xL=pa.getX(j)/lunit
             yL=pa.getY(j)/lunit
             zL=pa.getZ(j)/lunit
             tau2=xL*xL+yL*yL+zL*zL 
             p=ParticleD(pa.getPx(j)/un, pa.getPy(j)/un, pa.getPz(j)/un, pa.getEnergy(j)/un)
             if (tau2>10000):
                 ll_particles.add(p)
             else:
                 particles.add(p)
      # print particles.size(), ll_particles.size() 
      # jet from not LLP
      ktjet.buildJets(particles);
      # ktjet.printJets();
      jets=ktjet.getJetsSorted()
      # jet from LLP
      ktjet.buildJets(ll_particles);
      # ktjet.printJets();
      jetsLL=ktjet.getJetsSorted()
      if (len(jets)>0 and len(jetsLL)>0):
        # print "Jet Delta phi=",jets[0].phi() - jetsLL[0].phi()
        # print "  good=",jets[0].perp()," phi=",jets[0].eta(), " phi=",jets[0].phi()
        # print "  from LLP=",jetsLL[0].perp()," phi=",jetsLL[0].eta(), " phi=",jetsLL[0].phi()
        h1.fill(jets[0].perp())
        h2.fill(jetsLL[0].perp())
        h3.fill( abs(jets[0].phi() - jetsLL[0].phi()) )

   stat = file.getStat()
   cross=stat.getCrossSectionAccumulated()
   erro=stat.getCrossSectionErrorAccumulated();
   file.close()

lumi=nev/cross;
print "Lumi=%.3e pb"%lumi
print "Total cross section (pb)=",cross

c1 = HPlot("=HepSim=",600,700,1,2)
c1.visible(True)
c1.setRange(0,2000,0,1000)
c1.setMarginLeft(90)
c1.setNameX("p_{T}(lead. jet) [GeV]")
c1.setNameY("Events")

c1.cd(1,1)
c1.draw(h1)
c1.draw(h2)

# h1.toTable() 

l1=HLabel("σ=%.3e pb"%cross, 0.4, 0.75, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 15))
c1.add(l1)

l2=HLabel("=HepSim=",0.7,0.85, "NDC")
l2.setColor(Color.gray)
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l2)

c1.cd(1,2)
c1.setRange(0,5,0,500)
c1.setMarginLeft(90)
c1.setNameX("| φ_{j1} - φ_{j2}|")
c1.setNameY("Events")
c1.draw(h3)


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