Here is an example showing how to generate events with ISAJET and PGS.

#!/usr/bin/env python
#------------------------------------------------------------------------------
# File: isajet.py
# Description: Example showing how to generate ISAJET events using the
#              pgs (Pretty Good Simulator) Python module
# Created April 2001, Harrison B. Prosper
#------------------------------------------------------------------------------
# IMPORTS
#---------------------------------------------------------------
import os, sys
import pgs
import hntuple
from math    import *          # contains math functions
from physics import *          # physics utilities

#---------------------------------------------------------------    
# CONSTANTS
#---------------------------------------------------------------    
NEVENT = 2000                  # Number of events to generate
NPASS  = 200                   # Number of events that should pass cuts
NPSTEP = 20                    # "I'm alive!" print out frequency

NPHICAL  = 64
NETACAL  = 84
CONESIZE = 0.7

MAXJET   = 4
MAXELE   = 4
MAXPHO   = 4
MAXMUO   = 4

OPTION = "ISAJET"
ENERGY = 2000.0
FILE   = "so10"
CARDS  = FILE + ".cmd"         # Command file
TREE   = FILE + ".tree"        # File for tree listing
TABLE  = FILE + ".table"       # File for table listing
NTUP   = FILE + ".hbook"       # Ntuple output file


#---------------------------------------------------------------    
# FUNCTIONS
#---------------------------------------------------------------    
def flavor(x,f):
    return abs(f-x) < 0.5


#---------------------------------------------------------------    
# Create empty ntuple
#
# The object "tag" is a Python list
#
# IMPORTANT: The names defined here must be exactly the same as the
#            variables used below. This list of tags is used by the
#            ntuple capture method to store the correct variables
#            in the ntuple.
#
#---------------------------------------------------------------    
tag = ["NJET", "NPHO", "NELE", "NMUO", "HTJ", "HT", "MET","METPHI"]

# Jet variables

for i in xrange(1,MAXJET+1):
    s = "ETJ%d" % i; tag.append(s)
    
for i in xrange(1,MAXJET+1):
    s = "PHIJ%d" % i; tag.append(s)

for i in xrange(1,MAXJET+1):
    s = "ETAJ%d" % i; tag.append(s)

for i in xrange(1,MAXJET+1):
    s = "FLAVJ%d" % i; tag.append(s)

# Photon variables

for i in xrange(1,MAXPHO+1):
    s = "ETP%d" % i; tag.append(s)
    
for i in xrange(1,MAXPHO+1):
    s = "PHIP%d" % i; tag.append(s)

for i in xrange(1,MAXPHO+1):
    s = "ETAP%d" % i; tag.append(s)

# Electron variables

for i in xrange(1,MAXELE+1):
    s = "ETE%d" % i; tag.append(s)
    
for i in xrange(1,MAXELE+1):
    s = "PHIE%d" % i; tag.append(s)

for i in xrange(1,MAXELE+1):
    s = "ETAE%d" % i; tag.append(s)

# Muon variables

for i in xrange(1,MAXMUO+1):
    s = "ETM%d" % i; tag.append(s)
    
for i in xrange(1,MAXMUO+1):
    s = "PHIM%d" % i; tag.append(s)

for i in xrange(1,MAXMUO+1):
    s = "ETAM%d" % i; tag.append(s)


# Create an output stream (basically, a random access file) for the ntuple

out = hntuple.onstream(NTUP)

# Create an empty ntuple that is attached to
# the stream that we have just opened

nt  = out.ntuple("1","SO(10)",tag)
    
#---------------------------------------------------------------    
# MAIN PROGRAM
#---------------------------------------------------------------    

# SETUP

pgs.setOption(OPTION,CARDS)
pgs.setNumberEvents(NEVENT)
pgs.setConeSize(CONESIZE)
pgs.setNphiCal(NPHICAL)
pgs.setNetaCal(NETACAL)
pgs.setBeams('p','pbar',ENERGY)

# Create printer object that will print events to file in a table format

printer = pgs.printer(TABLE,NPASS)

# Create a pretty printer to print events to file in a tree format

ppfile   = pgs.prettyPrinter(TREE)
ppscreen = pgs.prettyPrinter()

# Initialize PGS. NOTE: This must be called after the above setup

pgs.initialize()

# Create empty vectors to receive reconstructed objects

ve    = pgs.vobject()
vm    = pgs.vobject()
vj    = pgs.vobject()
vp    = pgs.vobject()


# --------- START OF LOOP
    
event1 = 0
event2 = 0

mathError  = []
otherError = []

startTime = time()

while pgs.nextEvent():           # <--------- GENERATE EVENTS

    event1 = event1 + 1
        
    # "I'm alive" print out

    if (event1 % NPSTEP) == 0:
        
        # Get particles at root of trees

        p0 = pgs.particle(0)     # First particle in HEPEVT
        p1 = pgs.particle(1)     # Second particle in HEPEVT
        
        ppfile.write(p0,event1);   ppfile.write(p1)   # Tree printout (to file)
        ppscreen.write(p0,event1); ppscreen.write(p1) # To screen

        # Estimate and display amount of time left
        
        ImAlive(FILE, event1, event2, NPASS, startTime)

        

    # Simulate detector and reconstruct event using PGS
    #---------------------------------------------------------
    
    pgs.reconstruct()
       
    # Get reconstructed jets 
    
    pgs.jets(vj)
    
    # Get reconstructed electrons

    pgs.electrons(ve)

    # Get reconstructed muons

    pgs.muons(vm)

   # Get reconstructed muons

    pgs.photons(vp)

    # Get missing ET object and define some convenient symbols
        
    miss = pgs.misset()


    # NOTE: We MUST use the exactly the same names as defined in
    #       the "tag" sequence above. Python is CaSe SeNsiTivE!

     
    # --------- Apply minimal cuts
    #-----------------------------

    NJET  = len(vj)
    if NJET       < 1    : continue

    ETJ1  = vj[0].et()
    if ETJ1       < 10.0 : continue
    
    ETAJ1 = vj[0].eta()
    if abs(ETAJ1) >  2   : continue
    
    MET   = miss.et()
    METPHI= miss.phi()
    
    NELE  = len(ve)
    NMUO  = len(vm)
    NPHO  = len(vp)

    
    # Clear variables
        
    for field in tag:
        exec("%s = 0.0" % field)
        
    # Set jet variables

    for i in xrange(1,min(NJET,MAXJET)+1):
        exec("ETJ%d    = vj[%d-1].et()"  % (i,i))
        exec("PHIJ%d   = vj[%d-1].phi()" % (i,i))
        exec("ETAJ%d   = vj[%d-1].eta()" % (i,i))
        exec("FLAVJ%d  = vj[%d-1].item(pgs.rFLAVOR)"   % (i,i))

    # Set electron variables
    
    for i in xrange(1,min(NELE,MAXELE)+1):
        exec("ETE%d    = ve[%d-1].et()"  % (i,i))
        exec("PHIE%d   = ve[%d-1].phi()" % (i,i))
        exec("ETAE%d   = ve[%d-1].eta()" % (i,i))

    # Set muon variables
    
    for i in xrange(1,min(NMUO,MAXMUO)+1):
        exec("ETM%d    = vm[%d-1].et()"  % (i,i))
        exec("PHIM%d   = vm[%d-1].phi()" % (i,i))
        exec("ETAM%d   = vm[%d-1].eta()" % (i,i))

    # Set photon variables
    
    for i in xrange(1,min(NPHO,MAXPHO)+1):
        exec("ETP%d    = vp[%d-1].et()"  % (i,i))
        exec("PHIP%d   = vp[%d-1].phi()" % (i,i))
        exec("ETAP%d   = vp[%d-1].eta()" % (i,i))


    #---------------------------------------
    # Compute some other stuff
    #---------------------------------------

    HTJ = 0
    for jet in vj:
        HTJ = HTJ + jet.et()

    HT = HTJ
    for e in ve:
        HT  = HTJ + e.et()


    #---------------------------------------
    # Save variables
    #---------------------------------------

    nt.capture()


    event2 = event2 + 1
    if event2 >= NPASS : break
    
    # --------- END OF LOOP


ppfile.close()
printer.close()
out.close()

# Write summary

ImDone(FILE, event1, event2, startTime)