% 04/08/87 708051606 MEMBER NAME INFOL52 (LUND) M TEX C######################################################################C C C C L E P T O C C C C version 5.2, July 31, 1987 C C C C The Lund Monte Carlo for Deep Inelastic Lepton-Nucleon Scattering C C C C Author: Gunnar Ingelman, DESY theory group, C C Notkestrasse 85, D-2000 Hamburg 52, FRG C C EARN/BITNET: T00ING @ DHHDESY3 C C Contribution on parton cascades from C C M. Bengtsson, T. Sjoestrand, Lund University C C C C Documentation: short manual in T00ING.LUND(INFOL52) C C Source code in T00ING.LUND(LEPTO52) C C Please report any errors or suggestions for impromevents. C C C C######################################################################C 1. Introduction. =============== The current program, LEPTO 5.2, simulates complete deep inelastic lepton-nucleon scattering events based on the standard electroweak theory, including some QCD corrections, for the hard parton level interaction and the Lund string model {1} for the subsequent soft hadronization process. It is a further development of earlier versions {2,3}, but not backwards compatible with those. 2. Physics input. ================ Electron, muon, neutrino and their antiparticles are valid as incoming leptons. The target nucleon can either be a proton (default) or a 'nucleus' in the sense of a simple sum of nucleons such that, for each event a proton or a neutron is chosen for the interaction based on their probabilities in terms of summed parton structure functions. Note however, that no other nuclear effects are taken into account and that the remaining part of the 'nucleus' is neglected once a target nucleon has been selected. The momenta of incoming particles can be specified in any frame and options exist for transforming the final event into a few other standard frames. All electroweak interactions can be simulated, i.e. pure photon exchange, pure Z exchange, photon-Z exchange with interference and W+- exchange. The basic kinematical variables, like x and y or x and Q**2, are chosen according to the electroweak cross section formulae as given in {4} or, optionally, they can be provided by the user. By specifying cuts, a specific region in phase space can be defined and events generated in that region only. It should be noted that the cross section formulae used for this purpose are based on the pure quark-parton model with electroweak matrix elements for the hard scattering folded with QCD evolved quark structure functions. Thus, the longitudinal structure function and similar QCD corrections are not included in this version, in particular the contribution to F2 from heavy flavour production via boson-gluon fusion is only included in an approximate way through the possibility of a heavy quark structure function parametrization. The lepton polarization dependence is, however, included. Hence, the polarization of the incoming lepton beam can be specified and, for each event, a certain helicity state is chosen based on its cross section weighted with the probability for that helicity state. The total cross section for the specified process can be obtained in two ways. The first is by performing a 2-dimensional numerical integration in the initialization (see LST(10)) using an adaptive Gaussian method. The required accuracy can be specified, by PARL(15), and the result is printed and stored in PARL(23). Thus, if only the cross-section in a given kinematical region is needed, this can be obtained without the more time consuming event simulation. The second method is based on the Monte Carlo simulation of events which samples the differential cross section. The resulting estimate, stored in PARL(24), is updated with each generated event and its accuracy depends on the number of generated events as 1/sqrt(N). Note that, with both methods, the given cross section corresponds to the region allowed by the specified cuts only and can therefore be directly used to normalize the generated event sample. The first order QCD processes of gluon radiation from the struck quark and boson-gluon fusion into a quark-antiquark pair are included according to their proper matrix elements {5}. The latter, however, in an approximation where the quark masses are not explicitly included in the matrix element, but only taken into account through an overall threshold factor. To decide on an event by event basis whether to generate one of these event types (qg- or qq-event) rather than the leading order process (q-event), the probabilities for each event type must be available as a function of the kinematical variables, here chosen as x and W. This involves evaluation of 3-dimensional integrals over the extra degrees of freedom (like energy, polar and azimuthal angle of the quark) available in the three-body final state of the QCD processes compared to the two-body final state (quark and lepton) in the leading process. The integration over azimuthal angle is trivial and, using the variables zq and xp {5}, the zq-variable can be integrated analytically leaving only the xp-integration, which involves structure functions, to be performed numerically. To save time this is not done for each event but rather only once in the initialization using an adaptive Gaussian method and the corresponding qg- and qq-event probabilities then stored (in common LGRID) on a 'grid' in the x-W plane such that the probability for any x,W are obtained by interpolation when generating events. The information on the 'grid' can be saved on a file for use in later jobs, but due to the reduction to 1-dimensional numerical integration in the present version there is no strong need to do so, since setting up the grid only takes a few cpu seconds in the initialization. As an alternative to these exact, but fixed order perturbative approach, QCD parton cascade evolution, to effectively simulate higher order effects, is also included, see LST(8). Algorithms for the parton shower in initial and final state used in e+e- and hadron collider physics have been specially adopted for deep inelastic scattering, see {12} for details and some essential results. This approach involves approximations, essentially leading log, and is therefore not expected to give quite correct results for hard emissions giving clearly separable jets. On the other hand, the effective inclusion of higher order effects is expected to give a better description of multijet phenomena as well as the internal structure of a jet (softening and widening from multiple gluon emission). Unfortunately, there is not yet a prescription of how to take the first order exact QCD matrix elements properly into account in the parton cascade approach applied to deep inelastic scattering. Therefore, one has to make the choice of whether to follow the road of exact, but lowest order, matrix elements or that of parton showers. At fixed target energies, the two will not differ much in most observables, but at ep collider energies significant differences appear and the cascade approach will usually be preferable since it is expected to better simulate the complexity of the events. The target remnant, left over when a quark is struck, is in the simplest case a diquark. If, however, the struck quark was a sea quark this remnant will also contain the antiquark of the struck quark and thus be a more complicated 4-quark object which must be dealt with in order to conserve the flavour properties. If this remaining partner sea quark is an anti-up or anti-down, it is simply considered to annihilate a corresponding valence quark and a diquark system is again obtained. If it is an anti-strange quark, it will form a strange meson together with a valence quark chosen at random leaving a diquark that forms a string system with the struck s-quark. If an antiquark interacted, the left over quark partner will form a baryon together with a valence diquark leaving a randomly chosen valence quark that forms a string system with the kicked out antiquark. The fragmentation of a diquark is covered by the Lund model, see e.g. {6} for a study in DIS, but the formation of an 'extra' hadron from more complicated target remnant systems is not and must be modelled. This splitting of a target remnant into a hadron and a parton is made by giving them equal and opposite pt, PARL(14), like in the string breaking process and sharing the available longitudinal momentum based on counting rules. The partner sea quark is here assumed to have no essential dynamic property and hence the three valence quarks is split such that the quark takes an energy-momentum fraction, x, given by dP/dx = 2*(1-x) which implies that a valence quark on the average takes one third of the available energy-momentum. Also in the case of boson-gluon fusion a non-trivial target remnant consisting of the three valence quarks in a colour octet state occurs. This is treated in a similar way by splitting it into a quark and a diquark that forms separate colour triplet strings together with the produced antiquark and quark respectively. For more details on the model for baryon production we refer to {6}, for the model in general to {7} and for QCD effects to {7,8,12}. 3. Description of program components. ==================================== The program is written completely in FORTRAN 77 and is a further development of earlier versions of LEPTO {2,3}, but is not backwards compatible due to changed common blocks and subroutine calls. 3.1 Subroutines and functions. ----------------------------- Here we only describe the routines that the user will certainly call and a few that may be be useful, but leave the description of all internal routines for the coming full documentation. SUBROUTINE LINIT(LFILE,LEPIN,PLX,PLY,PLZ,PPX,PPY,PPZ,INTER) Purpose: to initialize the generation procedure. This includes finding maximum of that part of the differential cross section which is used for the rejection procedure, obtaining the information (stored in LGRID) for the QCD event probabilities and (optionally, LST(10)) calculating the cross section using an adaptive Gaussian method. Argumets: LFILE : logical file number for QCD weights. For positive values they are read from an existing file, whereas for negative values they are first calculated (if LST(8)=1) and then stored on file number -LFILE. If 0 no weights file is used, the weights are calculated if LST(8)=1 but not saved. Note: since it only takes some cpu seconds to calculate these weights, there is no strong reason to store them. LEPIN : is type of lepton (Lund code). PLX,PLY,PLZ : is the momentum vector for the lepton (in GeV/c). PPX,PPY,PPZ : is the momentum vector for the nucleon. INTER is the type of interaction to be simulated: = 1 : electromagnetic (EM) = 2 : weak charged current (CC) = 3 : weak neutral current = 4 : neutral current (NC), i.e. gamma+Z with interference SUBROUTINE LEPTO Purpose: to administer the generation of one event of the kind specified in the last LINIT call. SUBROUTINE LPRWTS(NSTEP) Purpose: can be called anytime after a call to LINIT to print the probabilities for q-, qg- and qq-event types using the present QCD weights in common LGRID. Only the values on each NSTEP point in the x,W grid is printed, i.e. NSTEP=1 prints the values on all grid points. SUBROUTINE LFRAME(IFRAME,IPHI) Purpose: to transform the event between different frames. Arguments: LFRAME: specifies the wanted frame (as for LST(5)): = 1 : hadronic CM frame, z-axis along boson. = 2 : lepton-nucleon CM frame, z-axis along lepton. = 3 : lab system as defined by the user. = 4 : as 3 but z-axis along exchanged boson. IPHI: specifies whether a random rotation in azimuthal angle phi of the scattering plane is to be made or not. = 0 : no rotation; if it has been done earlier such a rotation is undone such that the x-z plane is defined by the scattered lepton. = 1 : random rotation. Remaks: The phi-rotation is always performed in the lepton-nucleon CMS frame; if IFRAME = 1 it is never made. The present frame is stored in LST(28), LST(29) and is updated by a call to LFRAME. SUBROUTINE LWWB(ENU) Purpose: to give energy (ENU) of a (anti-)neutrino chosen from a simple parametrization (defined by DATA statement within the routine) of a wide band beam. The energy is actually chosen from the beam spectrum weighted with the energy since the cross section rises linearly with energy; an effect not taken into account in LEPTO. Note: this routine is only meant as an example and should rather be replaced by a more realistic beam energy distribution. SUBROUTINE PYSTFU(KF,X,Q2,XPQ) Purpose: to give values of parton structure functions. (This routine is obtained from PYTHIA 4.8 {9}). Arguments: KF: particle flavour code; 41=proton, 42=neutron, 17=pi+, negative values for the corresponding antiparticles. X: momentum carried by parton. Q2: momentum transfer. XPQ: array (-6:6) that on return contains the values of the parton structure functions. Index: 0=gluon, 1=u, 2=d, 3=s, 4=c, 5=b 6=t and negative values for antiquarks. Remark: The output is xq(x,Q2) etc. The parametrisation is chosen by LST(15) and max flavour can be restricted by LST(12). SUBROUTINE PYTIME(TIME) Purpose: to get the elapsed time by the call to some machine dependent routine. Default is the TIMEX routine in the CERN library. Remark: since this information is not essential, TIME can simply be set to zero if a suitable routine is not available. 3.2 Common blocks. ----------------- Most of the communication between the user and the program is via the switches and parameters in the common blocks. Essentially the user need only be concerned with common LEPTOU since all others are only for internal use. In block data LEPTOD all parameters are given sensible default values, indicated by (D=...) below. The generated event is stored in common block LUJETS, described in {10}, which the user must be familiar with. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U Purpose: contains basic switches and parameters which regulates the performance of the program. Also kinematical variables and parameters for applying cuts on them. Note that LST(2)-LST(19) and PARL(1)-PARL(19) are input switches and parameters, whereas LST(20)-LST(40) and PARL(20)-PARL(30) are output flags and parameters. Parameters: CUT(1), CUT(2): (D=.001,1.) lower and upper limit of Bjorken-x variable. CUT(3), CUT(4): (D=0.00,1.) lower and upper limit of y variable. CUT(5), CUT(6): (D=4.,1.E8) lower and upper limit of Q**2 (GeV**2). CUT(7), CUT(8): (D=5.,1.E8) lower and upper limit of W**2 (GeV**2). CUT(9), CUT(10): (D=1.,1.E8) lower and upper limit of nu variable (GeV). CUT(11), CUT(12): (D=1.,1.E8) lower and upper limit of scattered lepton energy (GeV, in frame defined by LINIT call). CUT(13), CUT(14): (D=0.,3.1416) lower and upper limit of scattered lepton angle (rad, in frame defined by LINIT call). Remarks: The cuts are applied already when chosing kinematic variables, before evaluating cross section formulae and structure functions, and will therefore be more efficient than applying cuts later in the users program. The cross section estimates provided apply to the region allowed by the cuts. LST(1) : interaction type as given by INTER in LINIT call. LST(2) : (D=1) choice of simulation and applying cuts. =1 : choose x, y from differential cross section and apply cuts defined by CUT. =2 : x,y supplied by user via LEPTOU, apply cuts. =3 : x,y supplied by user via LEPTOU, no cuts. LST(3) : (D=1) complete event simulation or only kinematics choice =0 : only Monte Carlo integration of cross section =1 : full event simulation LST(4) : (D=1) final lepton included in event record or not. =0 : deactivates the scattered lepton in the event record by setting K(4,1)=40000, LUEDIT will then deleted it and leave the hadronic part of the event. =1 : scattered lepton included as final state particle. Note: scattered lepton is in line 4 of the event record, unless parton cascades are used in which case it is in line 9 (both assuming no LUEDIT call made). LST(5) : (D=3) choice of frame for event. =1 : hadronic CM frame, z-axis along boson =2 : lepton-nucleon CM frame, z-axis along lepton =3 : lab system as defined by the user =4 : as 3 but z-axis along exchanged boson LST(6) : (D=1) choice of performing random rotation for azimuthal angle of scattering plane defined by scattered lepton. =0 : no phi-rotation, scattering plane is x-z plane. =1 : random phi-rotation, scattering plane is x-z plane. Note that this rotation is made in the lepton-nucleon cms frame LST(7) : (D=1) choice of performing fragmentation =0 : skip fragmentation, only parton level event generated. (Hadronization can be done later by calling LUEXEC.) =1 : fragmentation performed. LST(8) : (D=1) simulation of QCD events =0 : QCD events switched off =1 : 1:st order QCD processes (gluon radiation and boson- gluon fusion) according to exact matrix elements. =2 : QCD parton cascade evolution of initial and final quark =3 : QCD parton cascade evolution of initial quark only. =4 : QCD parton cascade evolution of final quark only. =5 : QCD swithced off, but target treatment exactly as in parton cascade case. LST(9) : (D=4) max no. of flavours in expression for alpha-strong LST(10): (D=1) choice of making adaptive Gaussian integration of cross section in LINIT call. =0 : excluded =1 : included Note, that the integration region is only that allowed by the cuts, required accuracy specified by PARL(15) and the result is stored in PARL(23). LST(11): (D=4) lowest quark flavour for which a threshold factor is applied in boson-gluon fusion into q+qbar. Note: the heavy quark mass is not explicitly included in the matrix element, only simple threshold factor is applied. LST(12): (D=4) max flavour in sea structure functions. LST(13): (D=4) heaviest quark flavour that can be produced in boson-gluon fusion process LST(14): (D=1) baryon production from target remnant =0 : excluded, i.e. target remnant typically approximated by an antiquark instead of a diquark. =1 : included. LST(15): (D=1) choice of proton structure function parametrisations = 0: simple scaling functions = 1: EHLQ set 1 {13} = 2: EHLQ set 2 {13} = 3: Duke-Owens set 1 {14} = 4: Duke-Owens set 2 {14} = 5: Gluck-Hoffman-Reya {15} LST(16): (D=0) choice of variable beam energy of incoming particles =0: energy not allowed to vary from event to event. =1: energy allowed to vary, should be given in P(i,j) with i=1,2 for lepton, nucleon and j=1..5 for px,py,pz,E,mass. Note: this option is not well tested yet and is not operationally with first order QCD switched on. LST(17): (D=2) max virtuality in parton cascades given by =1: Q**2 =2: W**2 LST(18): (D=0) concerning W and Z masses =0: masses calculated from standard model using sin(theta), alpha-em, GF and radiative corrections; i.e. PARL(5), PARL(16)--PARL(18) =1: Z, W masses taken from PMAS(2), PMAS(3) in common LUDAT2 and can be varied independently of other variables. (Note: default masses not consistent with previous option.) LST(19): (D=3) choice of grid for 1st order QCD weights =0: user defined grid, read in free format by subroutine LWEITS =1: grid suitable for fixed target, beam energy < 300 GeV =2: grid suitable for fixed target, beam energy <1000 GeV =3: grid suitable for HERA =4: grid suitable for LEP+LHC LST(20): Internal flag, 0/1 = simulation integration. LST(21): error flag, should normally be zero to indicate that an event is properly generated. Nonzero values indicates an error in which case no event was generated. This may occur in case of variable beam energy and user supplied kinematical variables. LST(22): specifies chosen target nucleon in current event =1: proton =2: neutron LST(23): specifies process simulated =1: EM interaction =2: W emitted from lepton =3: W emitted from antilepton =4: Z0 emitted from lepton =5: Z0 emitted from antilepton =6: EM + NC + interference LST(24): specifies current event wrt 1st order QCD matrix element =1: q-jet event, i.e. leading order process. =2: qg-jet event, i.e. gluon radiation event. =3: qq-jet event, i.e. boson-gluon into q+qbar event. Note: not set when using parton showers. LST(25): specifies flavour of struck quark in current event; Lund IFL code, 1=u, 2=d etc. LST(26): gives entry line (in event record) of outgoing struck quark Note: not set when using parton showers. LST(27): flag to signal extra hadron from target remnant =0: no extra hadron formed. =1: extra hadron formed due to non-trivial target remnant, e.g. remaining 4-quark state when sea antiquark was struck. Note: not set when using parton showers. LST(28): flag to tell what frame the event is given in =1: hadronic CM frame, z-axis along boson =2: lepton-nucleon CM frame, z-axis along lepton =3: lab system as defined by the user =4: as 3 but z-axis along exchanged boson LST(29): flag to specify whether a rotation in azimuthal angle, phi, of the scattering plane, defined by scattered lepton, has been performed or not =0: no phi-rotation. =1: random phi-rotation made. LST(30): gives helicity state of interacting lepton in current event, =-1: left handed =+1: right handed LST(31): specifies in what variables the simulation is performed =0: x and Q**2 =1: x and y LST(32)-LST(40): unused at present PARL(1) : (D=1.) number of nucleons in target nucleus PARL(2) : (D=1.) number of protons in target nucleus PARL(3) : (D=0.44 GeV) width of gaussian primordial transverse momentum distribution. PARL(4) : (D=0.75) probability that a ud-diquark in target remnant has spin and isospin equal zero, i.e. I=S=0. PARL(5) : (D=0.226) sin**2(theta-weinberg) PARL(6) : (D=0.) polarization (helicity) of incomimg lepton defined such that -1(+1) corresponds to full left-(right)handed beam. PARL(7) : (D=4. GeV**2) minimum Q**2 in alpha-strong below which alpha-strong is kept constant. PARL(8),PARL(9): (D=0.0025,2. GeV) minimum value of y=m_ij**2/W**2 and m_ij resp, where m_ij is the invariant mass of any parton pair (including target remnant), used as cutoffs against divergences in 1st order QCD matrix elements. Note: PARL(8) is the dominating cut at large energy (W) whereas PARL(9) dominates at low energy. These are starting values when integrating first order QCD matrix elements, but the effective cut used (PARL(27)) is automatically increased in order that the QCD probabilities do not exceed unity for any (x,W)-value. PARL(10): (D=0.3 GeV) lambda-QCD in alpha-s for 1st order QCD. PARL(11): (D=0.01) required relative accuracy when integrating QCD 1st order matrix elements. PARL(12), PARL(13): (D=4.,0.1) internal parameters for matrix element integration. PARL(14) (D=0.44 GeV) width of Gaussian pt when non-trivial target remnant is split into a particle and a jet. PARL(15) (D=0.01) required relative accuracy in adaptive Gaussian integration of cross section performed (optionally) in LINIT. PARL(16) (D=7.2973E-03) alpha-em PARL(17) (D=1.16637E-05 GeV**-2) weak Fermi constant PARL(18) (D=0.0711) delta-R for radiative correction for calculating W, Z boson masses. PARL(19)-PARL(20): unused at present PARL(21): 2Pk, i.e. twice the scalar product of incoming proton and lepton 4-vectors, equals s (inv. mass squared) if masses are neglected. PARL(22): 2Pq, i.e. twice the proton-boson scalar product. PARL(23): estimate of cross section in pb calculated in LINIT, see LST(10)=1, using an adaptive Gaussian integration routine. Only the kinematic region allowed by the cuts is considered. PARL(24): Monte Carlo estimate of cross-section (pb) associated with the generated event sample (taking cuts into account). Reset to at initialization and updated with each generated event. PARL(25): value of alpha-strong in current event. PARL(26): value lambda-QCD used in last structure function call. PARL(27): present value of y-cut for 1:st order QCD, given by maximum of PARL(8) and PARL(9), but also modified internally to obtain QCD event probabilities not exceeding unity. PARL(28)-PARL(30): unused at present X: Bjorken-x, x=Q**2/(2Pq) Y: normal y-variable, y=Pq/Pk W2: mass-squared of the hadronic system, W**2=(P+q)**2=Q**2 (1-x)/x +m**2 Q2: momentum tranfer squared, Q**2=-q**2 U: energy transfer, nu = Pq/m Here P, k, q are the four-vectors of the incoming nucleon, lepton and exchanged boson, respectively, and m is the nucleon mass. The following common blocks are used internally: COMMON /LFLMIX/ CABIBO(4,4) Purpose: Kobayashi-Maskawa matrix elements squared for flavour mixing. First index correspond to u, c, t, h flavours and second index to d, s, b, l quarks. Default values, from {11}, are CABIBO/.95,.05,2*0.,.05,.948,.002,2*0.,.002,.998,4*0.,1./ COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ COMMON /LINTEG/ NTOT,NPASS Purpose: internal parameters, charges and electroweak couplings, effective limits on kinematical variables. COMMON /LBOOST/ DBETA(2,3),DTHETA(2),DPHI(2),DE,DPZ,DPT,PB(5),PHIR Purpose: rotation and boost parameters. COMMON /LGRID/ NXX,NWW,XX(15),WW(15),PQG(15,15,3),PQQB(15,15,2), &QGMAX(15,15,3),QQBMAX(15,15,2),YCUT(15,15),XTOT(15,15),NP Purpose: information on grid for 1st order QCD probabilities. COMMON /LOPTIM/ OPTX(4),OPTY(4),OPTQ2(4),COMFAC Purpose: parameters to optimize speed of program. COMMON /PYMINU/ XKIN(4),UKIN(4),WKIN(4),AIN(4),BIN(4), &MAXFIN,RELUP,RELERR,RELER2,FCNMAX COMMON /PYMINC/ NAMKIN(4),NAM(30) Purpose: parameters for finding maximum of differential cross section. COMMON /PYPARA/ IPY(80),PYPAR(80),PYVAR(80) Purpose: parameters for parton cascade routines. This common block, which is from PYTHIA 4.8, is describe in detail in {9}. Only those parameters actually used with LEPTO in parton cascade mode are commented here (default values as in PYTHIA if not given): IPY(8) set to LST(12) in LINIT. IPY(13) set to zero in LINIT if LST(8)=3 or 5 IPY(14) set to zero in LINIT if LST(8)=4 or 5 IPY(15), IPY(16), IPY(17), IPY(34), IPY(40), IPY(41), IPY(42), IPY(47), IPY(48) are used. IPY(77) set to LST(17) in LINIT. PYPAR(7) set equal PARL(3) in LINIT. PYPAR(8), PYPAR(11) -- PYPAR(16) are used. PYPAR(21) (D=0.25 GeV) Lambda-QCD in initial state parton shower. PYPAR(22) (D=1. GeV**2) cutoff for initial state parton shower. PYPAR(23), PYPAR(24) are used. PYPAR(25), PYPAR(26) (D=2*1.) multiplies the chosen Q**2 scale to give maximum virtuality for final and initial state showers, respectively. PYPAR(27) (D=1.) multiplies virtuality Q**2 for alpha-s and structure functions in initial state showers. PYVAR(1) --- PYVAR(5) are used. 4. Examples on how to use the program ===================================== The Monte Carlo is a slave system, i.e. the user supplies his own main program to call the LEPTO routines. The program has to be loaded with the fragmentation routines JETSET 6.3 {10}. Two external functions are required. RLU to provide a uniformly distributed random number between zero and one; a dummy interface routine is provided in JETSET with examples for some computers. The ordinary gamma function is also called, GAMMA(X), but is often available in FORTRAN 77. The optional routine PYTIME to give the elapsed cpu time is also called and uses by default the TIMEX routine in the CERN library, but can be replaced or made into a dummy routine since this information is not essential. Below follows a few complete sample jobs to illustrate the use of the program. Example 1: Storing QCD weights on file. -------------------------------------- // JOB CLASS=E,TIME=(0,15),MSGLEVEL=(2,0) //*MAIN ORG=EXT,RELPRI=MED // EXEC JFORTCG,LLB1='R01UTL.CERN.PACKLIB' //C.SYSIN DD * C...Calculate QCD-weights for electron-proton neutral current C...interaction at HERA energy and store them on file number LFILE. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U C...1:st order QCD included, logical file number to store weights on. LST(8)=1 LFILE=20 C...Skip cross section integration. LST(10)=0 CALL LINIT(-LFILE,7,0.,0.,30.,0.,0.,-820.,4) END //G.SYSLIN DD // DD DSN=T00ING.OBJECT(LEPTO52),DISP=SHR // DD DSN=T00ING.OBJECT(JETSET63),DISP=SHR //G.FT20F001 DD DSN=T00ING.HERAW.NC, // DISP=(NEW,CATLG),UNIT=FAST,SPACE=(TRK,(1,1),RLSE) Example 2: Reading QCD weights from file and simulating events. -------------------------------------------------------------- // JOB CLASS=E,TIME=(0,15),MSGLEVEL=(2,0) //*MAIN ORG=EXT,RELPRI=MED // EXEC JFORTCG,LLB1='R01UTL.CERN.PACKLIB' //C.SYSIN DD * C...Read QCD-weights from file number LFILE and simulate C...electron-proton neutral current events at HERA energy. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LUJETS/ N,K(2000,2),P(2000,5) DIMENSION JETYPE(0:3) DATA JETYPE/4*0/ C...Cuts on kinematical variables. READ(5,*) CUT C...1:st order QCD included, logical file number to read weights from. LST(8)=1 LFILE=20 CALL LINIT(LFILE,7,0.,0.,30.,0.,0.,-820.,4) CALL LPRWTS(1) DO 100 NE=1,1000 10 CALL TIMEL(TIML) IF(TIML.LT.2.) GOTO 110 CALL LEPTO IF(LST(21).NE.0) GOTO 10 C...# of events, total and q-, qg-, qq-events; list first of each. JETYPE(0)=JETYPE(0)+1 JETYPE(LST(24))=JETYPE(LST(24))+1 IF(JETYPE(LST(24)).EQ.1) CALL LULIST(11) 100 CONTINUE 110 CONTINUE WRITE(6,1000) JETYPE,(100.*JETYPE(I)/FLOAT(JETYPE(0)),I=0,3) 1000 FORMAT(//,' Events generated:',3X,'total',7X,'q',6X,'qg',6X,'qq', &/,17X,'#',4I8,/,17X,'%',4F8.1) END /* //G.SYSLIN DD // DD DSN=T00ING.OBJECT(LEPTO52),DISP=SHR // DD DSN=T00ING.OBJECT(JETSET63),DISP=SHR //G.SYSIN DD * 0.01,1.,0.,1.,1000.,1.E+06,5.,1.E+06,1.,1.E+06,1.,1.E+06,0.,3.14159 /* //G.FT20F001 DD DSN=T00ING.HERAW.NC,DISP=SHR Example 3: Simulating events and storing them on file. ----------------------------------------------------- // JOB CLASS=A,TIME=(1,00),MSGLEVEL=(2,0) //*MAIN ORG=EXT,RELPRI=MED // EXEC JFORTCG,LLB1='R01UTL.CERN.PACKLIB' //C.SYSIN DD * C...Generate HERA events (with LEPTO 5.2) and store them on file. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LBOOST/ DBETA(2,3),DTHETA(2),DPHI(2),DE,DPZ,DPT,PB(5),PHIR COMMON /LUJETS/ N,K(2000,2),P(2000,5) DOUBLE PRECISION DTHETA,DPHI,DBETA,DE,DPZ,DPT CHARACTER*80 COMMNT(10) DATA NEV/0/ C...Max # events to be generated. NEVENT=10000 C...Lepton type: 7->e-, -7->e+. LEPIN=7 C...Lepton and proton momentum (along z-axis, with sign). PLEPIN=30. PPIN=-820. C...Interaction: 1->EM, 2->CC, 3->Weak NC, 4->NC (EM+weak) INTER=4 C...Lepton polarisation = (left-right)/(left+right), i.e. 0-> no pol. C...-1->full left handed, +1->full right handed polarisation. PARL(6)=0. C...QCD events: 0->off, 1->1st order matrix elements, 2->parton cascades LST(8)=2 C...Logical file # for QCD weights (for read/write). LFILE=0 C...Logical file # for storing events. NEVFIL=30 C...Comment lines (10 lines with 80 characters) to specify event file. READ(5,1100) COMMNT C...Specify kinematical cuts. READ(5,*) CUT C...Set seed for RN random number generator; default is standard seed. CALL RNSAVE(ISEED) CALL RNSET(ISEED) WRITE(6,1100) COMMNT WRITE(6,*) ' # events requested: ',NEVENT WRITE(6,*) ' Lepton type: ',LEPIN WRITE(6,*) ' Lepton momentum: ',PLEPIN WRITE(6,*) ' Proton momentum: ',PPIN WRITE(6,*) ' Type of interaction: ',INTER WRITE(6,*) ' QCD switch, LST(8) = ',LST(8) WRITE(6,*) ' File # for weights: ',LFILE WRITE(6,*) ' File # for events : ',NEVFIL WRITE(6,*) ' Starting random number seed = ',ISEED CALL TIMEX(TIME1) CALL LINIT(LFILE,LEPIN,0.,0.,PLEPIN,0.,0.,PPIN,INTER) C...Print weights if read from file. IF(LFILE.GT.0) CALL LPRWTS(1) CALL TIMEX(TIME2) WRITE(NEVFIL) COMMNT,LEPIN,PLEPIN,PPIN,INTER,NEVENT,CUT DO 500 IEV=1,NEVENT 100 CALL TIMEL(TIML) IF(TIML.LT.2.) GOTO 600 CALL LEPTO IF(LST(21).NE.0) GOTO 100 NEV=NEV+1 IF(NEV.LE.2) CALL LULIST(11) CALL RNSAVE(ISEED) C...Save on file: kinematics, parameters, seed, boosts, event record. WRITE(NEVFIL) 1,N,NEV,X,Y,W2,Q2,U,LST,PARL,ISEED WRITE(NEVFIL) 2,DBETA,DTHETA,DPHI,DE,DPZ,DPT,PB,PHIR WRITE(NEVFIL) 3,((K(I,J),I=1,N),J=1,2),((P(I,J),I=1,N),J=1,5) 500 CONTINUE 600 CONTINUE CALL RNSAVE(ISEED) WRITE(NEVFIL) 0,N,NEV,X,Y,W2,Q2,U,LST,PARL,ISEED CALL TIMEX(TIME3) WRITE(6,*) WRITE(6,*) ' # events generated :',NEV WRITE(6,*) ' Final random number seed = ',ISEED WRITE(6,*) WRITE(6,*) ' Time for initialization = ',TIME2-TIME1,' seconds.' WRITE(6,*) ' Time for event generation = ',TIME3-TIME2,' seconds.' WRITE(6,*) ' Time needed for one event = ', &(TIME3-TIME2)/MAX(1,NEV),' seconds.' 1100 FORMAT(A80) END //G.SYSLIN DD // DD DSN=T00ING.OBJECT(LEPTO52),DISP=SHR // DD DSN=T00ING.OBJECT(JETSET63),DISP=SHR //G.SYSIN DD * comment line 1: LEPTO version 5.2 comment line 2: HERA NC events (gamma+Z exchange with interference) comment line 3: e-, pz = 30 GeV, no polarization comment line 4: proton, pz = -820 GeV comment line 5: Q**2 > 1000 GeV**2, x > 0.01, W**2 > 10. GeV**2 comment line 6: 1:st order QCD (qg and qqb events) switched on comment line 7: 10000 events requested (less may have been generated) comment line 8: events stored in the HERA frame (but boost parameters comment line 9: saved so transformations via LFRAME possible later) comment line 10: 0.01,1.,0.,1.,1000.,1.E+06,10.,1.E+06,1.,1.E+06,1.,1.E+06,0.,3.14159 //* QCD weights can be stored on or read from file # 20. //*G.FT20F001 DD DSN=T00ING.HERAW.NC, //* DISP=(NEW,CATLG),UNIT=FAST,SPACE=(TRK,(1,1),RLSE) //* DISP=SHR //* Generated events stored on file # 30. Approx. 80 tracks/1k events //G.FT30F001 DD DSN=T00ING.HERAEVT.NC1,DISP=(NEW,CATLG), // UNIT=FAST,SPACE=(TRK,(100,50),RLSE),DCB=R01DCB.VBS Example 4: Reading events from file and making analysis on them. --------------------------------------------------------------- // JOB CLASS=E,TIME=(0,15),MSGLEVEL=(2,0) //*MAIN ORG=EXT,RELPRI=MED // EXEC JFORTCLG,LRGN=1000K,LLB1='R01UTL.CERN.PACKLIB' //C.SYSIN DD * C...Read HERA events from file and analyse them; C...kinematics reconstruction using Jacquet-Blondel method. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LBOOST/ DBETA(2,3),DTHETA(2),DPHI(2),DE,DPZ,DPT,PB(5),PHIR COMMON /INIT/ PLEPIN,PPIN,INTER,THEMIN DOUBLE PRECISION DTHETA,DPHI,DBETA,DE,DPZ,DPT CHARACTER*80 COMMNT(10) DATA NEV/0/ C...Size of beam hole (radians) for Jacquet-Blondel analysis. THEMIN=0.07 C...Logical file # to read events from. NEVFIL=30 READ(NEVFIL) COMMNT,LEPIN,PLEPIN,PPIN,INTER,NEVENT,CUT WRITE(6,1100) COMMNT WRITE(6,*) ' Lepton type : ',LEPIN WRITE(6,*) ' Lepton momentum : ',PLEPIN WRITE(6,*) ' Proton momentum : ',PPIN WRITE(6,*) ' Interaction type : ',INTER WRITE(6,*) ' # events requested : ',NEVENT WRITE(6,*) ' Beam hole size : ',THEMIN WRITE(6,1000) CUT CALL BOOK DO 500 IEV=1,NEVENT CALL TIMEL(TIML) IF(TIML.LT.2.) GOTO 600 READ(NEVFIL,END=600,ERR=900) &IREC,N,NEVNT,X,Y,W2,Q2,U,LST,PARL,ISEED IF(IREC.EQ.0) GOTO 600 IF(IREC.NE.1) GOTO 500 READ(NEVFIL,END=600,ERR=910) &IREC,DBETA,DTHETA,DPHI,DE,DPZ,DPT,PB,PHIR IF(IREC.NE.2) GOTO 500 READ(NEVFIL,END=600,ERR=920) &IREC,((K(I,J),I=1,N),J=1,2),((P(I,J),I=1,N),J=1,5) IF(IREC.NE.3) GOTO 500 C...Event ready for analysis. NEV=NEV+1 IF(NEV.LE.1) CALL LULIST(11) CALL ANALYS(NEV) 500 CONTINUE 600 CONTINUE WRITE(6,*) WRITE(6,*) ' # events read: ',NEV WRITE(6,*) CALL FINALE(NEV) STOP 900 WRITE(6,*) ' Error in read, 1.' STOP 910 WRITE(6,*) ' Error in read, 2.' STOP 920 WRITE(6,*) ' Error in read, 3.' STOP 1000 FORMAT(' Cuts on kinematical variables; x, y, Q**2, W**2, nu,' &' E-lepout, theta-lepout',7(/,2F12.3)) 1100 FORMAT(A80) END C***************************************************************** SUBROUTINE BOOK C...Book histograms for HBOOK. COMMON / / HMEMOR(5000) CALL HLIMIT(5000) CALL HBOOK1(1,'1/Nev*dN/d(log(Q**2)) $',100,0.,5.) CALL HBOOK1(2,'1/Nev*dN/dx $',100,0.,1.) CALL HBOOK1(3,'1/Nev*dN/dy $',100,0.,1.) CALL HBOOK1(4,'(Q2R-Q2)/Q2 no beam hole$',100,-1.,1.) CALL HBOOK1(5,'(XR-X)/X no beam hole$',100,-1.,1.) CALL HBOOK1(6,'(YR-Y)/Y no beam hole$',100,-1.,1.) CALL HBOOK1(7,'(Q2R-Q2)/Q2 with beam hole$',100,-1.,1.) CALL HBOOK1(8,'(XR-X)/X with beam hole$',100,-1.,1.) CALL HBOOK1(9,'(YR-Y)/Y with beam hole$',100,-1.,1.) RETURN END C***************************************************************** SUBROUTINE ANALYS(NEV) C...User analysis routine, examplified by Jacquet-Blondel analysis. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /INIT/ PLEPIN,PPIN,INTER,THEMIN DATA PI/3.1415926/ CALL HF1(1,ALOG10(Q2),1.) CALL HF1(2,X,1.) CALL HF1(3,Y,1.) C...Keep only stable particles in event record, neutrinos skipped. CALL LUEDIT(2) IF(NEV.LE.1) CALL LULIST(11) C...Exclude scattered lepton in loops below. ILOW=2 IF(INTER.EQ.2) ILOW=1 SGN=SIGN(1.,PLEPIN) YR=0. PXSUM=0. PYSUM=0. NHADR=0 DO 100 I=ILOW,N C...Loop over particles, full angular coverage. NHADR=NHADR+1 YR=YR+P(I,4)+SGN*P(I,3) PXSUM=PXSUM+P(I,1) PYSUM=PYSUM+P(I,2) 100 CONTINUE IF(NHADR.EQ.0) RETURN YR=YR/(2.*ABS(PLEPIN)) Q2R=(PXSUM**2+PYSUM**2)/(1.-YR) XR=Q2R/(YR*4.*ABS(PLEPIN*PPIN)) CALL HF1(4,(Q2R-Q2)/Q2,1.) CALL HF1(5,(XR-X)/X,1.) CALL HF1(6,(YR-Y)/Y,1.) YR=0. PXSUM=0. PYSUM=0. NHADR=0 DO 200 I=ILOW,N C...Loop over particles, keep only those outside beampipe (THEMIN). THETA=PLU(I,13) IF(THETA.GT.PI/2.) THETA=PI-THETA IF(THETA.LT.THEMIN) GOTO 200 NHADR=NHADR+1 YR=YR+P(I,4)+SGN*P(I,3) PXSUM=PXSUM+P(I,1) PYSUM=PYSUM+P(I,2) 200 CONTINUE IF(NHADR.EQ.0) RETURN YR=YR/(2.*ABS(PLEPIN)) Q2R=(PXSUM**2+PYSUM**2)/(1.-YR) XR=Q2R/(YR*4.*ABS(PLEPIN*PPIN)) CALL HF1(7,(Q2R-Q2)/Q2,1.) CALL HF1(8,(XR-X)/X,1.) CALL HF1(9,(YR-Y)/Y,1.) RETURN END C***************************************************************** SUBROUTINE FINALE(NEV) C...Normalise and print histograms. CALL HOPERA(1,'+',1,1,100./5./MAX(1,NEV),0.) CALL HOPERA(2,'+',2,2,100./1./MAX(1,NEV),0.) CALL HOPERA(3,'+',3,3,100./1./MAX(1,NEV),0.) CALL HLOGAR(1) DO 100 I=4,9 CALL HOPERA(I,'+',I,I,100./2./MAX(1,NEV),0.) 100 CALL HLOGAR(I) CALL HPRINT(0) RETURN END //L.SYSLIN DD // DD DSN=T00ING.OBJECT(LEPTO52),DISP=SHR // DD DSN=T00ING.OBJECT(JETSET63),DISP=SHR //G.FT30F001 DD DSN=T00ING.HERAEVT.NC1,DISP=SHR Example 5: GEP-job cards. ------------------------ The above jobs can also be executed as GEP jobs at the DESY IBM with the following set of job cards: // JOB CLASS=A,TIME=(1,00),MSGLEVEL=(2,0) //*MAIN ORG=EXT,RELPRI=MED // EXEC GEP,LA=FORT77,COREGN=1200K, // VS=HBOOK9,OPT=3,FORT='GOSTMT,XL',REGION.GO=1500K, // DN='T00ING.GEP.TEST01',LIB0='SYS1.FT77LIB' C...FORTRAN 77 code as in above jobs. Note, however, that when using C...HBOOK the following two lines should not be included: C COMMON / / HMEMOR(5000) C CALL HLIMIT(5000) //LKED.SYSLIB DD // DD // DD // DD // DD // DD DSN=R01UTL.CERN.PACKLIB,DISP=SHR //LKED.PL DD DSN=T00ING.OBJ,DISP=SHR //LKED.SYSIN DD * INCLUDE PL(LEPTO52) INCLUDE PL(JETSET63) //GO.SYSIN DD * comment line 1: LEPTO version 5.2 comment line 2: HERA NC events (gamma+Z exchange with interference) comment line 3: e-, pz = 30 GeV, no polarization comment line 4: proton, pz = -820 GeV comment line 5: Q**2 > 1000 GeV**2, x > 0.01, W**2 > 10. GeV**2 comment line 6: 1:st order QCD (qg and qqb events) switched on comment line 7: 10000 events requested (less may have been generated) comment line 8: events stored in the HERA frame (but boost parameters comment line 9: saved so transformations via LFRAME possible later) comment line 10: 0.01,1.,0.,1.,1000.,1.E+06,10.,1.E+06,1.,1.E+06,1.,1.E+06,0.,3.14159 /* //* QCD weights can be written to or read from file # 20. //*GO.FT20F001 DD DSN=T00ING.HERAW.NC, //* DISP=(NEW,CATLG),UNIT=FAST,SPACE=(TRK,(5,5),RLSE) //* DISP=SHR //* Generated events stored on file # 30. Approx. 100 tracks/1k events //GO.FT30F001 DD DSN=T00ING.HERAEVT.NC,DISP=(NEW,CATLG), // UNIT=FAST,SPACE=(TRK,(100,50),RLSE),DCB=R01DCB.GEP References. ========== {1} B. Anderson, G. Gustafson, G. Ingelman, T. Sjostrand, Physics reports 97 (1983) 31 {2} G. Ingelman, T. Sjostrand, 'A Monte Carlo for leptoproduction', Lund preprint LU TP 80-12 {3} G. Ingelman, LEPTO version 4.3, CERN program pool long writeup program W5046, June 1986. {4} DESY HERA 80/01 ('the green book') {5} R.D. Peccei, R. Rueckl, Nucl. Phys. B162 (1980) 125 {6} B. Andersson, G. Gustafson, G. Ingelman, T. Sjostrand, Z. Phys. C13 (1982) 361 {7} B. Andersson, G. Gustafson, G. Ingelman, T. Sjostrand, Z. Phys. C9 (1981) 233 {8} G. Ingelman, B. Andersson, G. Gustafson, T. Sjostrand, Nucl. Phys. B206 (1982) 239 {9} H.U. Bengtsson, T. Sjostrand, 'The Lund Monte Carlo for hadronic processes', Lund preprint LU TP 87-3 {10} T. Sjostrand, 'The Lund Monte Carlo for jet fragmentation', JETSET version 6.2, LU TP 85-10, Comp. Phys. Commun. 39 (1986) 347 T. Sjostrand, M. Bengtsson, JETSET version 6.3 (updates), Lund preprint LU TP 86-22 {11} Review of particle properties, Phys. Lett. 170B (1986) 74 {12} M. Bengtsson, T. Sjostrand, LU TP 87-10 M. Bengtsson, G. Ingelman, T. Sjostrand, DESY preprint 87-??? {13} E. Eichten, I. Hinchliffe, K. Lane, C. Quigg, Rev. Mod. Phys. 56 (1984) 579, ibid. 58 (1986) 1047 {14} D.W.Duke, J.F.Owens, Phys. Rev. D30 (1984) 49 {15} M. Glueck, E. Hoffman, E. Reya, Z. Phys. C13 (1982) 119