Removed until updated.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 5 Nov 2008 10:03:35 +0000 (10:03 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 5 Nov 2008 10:03:35 +0000 (10:03 +0000)
463 files changed:
ISAJET/AliIsajetRndm.cxx [deleted file]
ISAJET/AliIsajetRndm.h [deleted file]
ISAJET/CMakeLists.txt [deleted file]
ISAJET/CMake_libisajet.txt [deleted file]
ISAJET/code/alqcd.F [deleted file]
ISAJET/code/amass.F [deleted file]
ISAJET/code/amgmw.F [deleted file]
ISAJET/code/charge.F [deleted file]
ISAJET/code/ctxc2i.F [deleted file]
ISAJET/code/ctxi2c.F [deleted file]
ISAJET/code/ctxin.F [deleted file]
ISAJET/code/ctxout.F [deleted file]
ISAJET/code/datime.F [deleted file]
ISAJET/code/dblpcm.F [deleted file]
ISAJET/code/dblvec.F [deleted file]
ISAJET/code/dboost.F [deleted file]
ISAJET/code/decay.F [deleted file]
ISAJET/code/decjet.F [deleted file]
ISAJET/code/decps1.F [deleted file]
ISAJET/code/decps2.F [deleted file]
ISAJET/code/decss3.F [deleted file]
ISAJET/code/dectau.F [deleted file]
ISAJET/code/decva.F [deleted file]
ISAJET/code/dhelas.F [deleted file]
ISAJET/code/dincgm.F [deleted file]
ISAJET/code/domssm.F [deleted file]
ISAJET/code/drllyn.F [deleted file]
ISAJET/code/ebeam.F [deleted file]
ISAJET/code/eebeg.F [deleted file]
ISAJET/code/eemax.F [deleted file]
ISAJET/code/elctrn.F [deleted file]
ISAJET/code/epf.F [deleted file]
ISAJET/code/estruc.F [deleted file]
ISAJET/code/evol01.F [deleted file]
ISAJET/code/evol02.F [deleted file]
ISAJET/code/evol03.F [deleted file]
ISAJET/code/evol05.F [deleted file]
ISAJET/code/evol06.F [deleted file]
ISAJET/code/evol07.F [deleted file]
ISAJET/code/evol11.F [deleted file]
ISAJET/code/evolms.F [deleted file]
ISAJET/code/evolve.F [deleted file]
ISAJET/code/fbrbm.F [deleted file]
ISAJET/code/flavor.F [deleted file]
ISAJET/code/fortop.F [deleted file]
ISAJET/code/frgjet.F [deleted file]
ISAJET/code/frgmnt.F [deleted file]
ISAJET/code/gamma.F [deleted file]
ISAJET/code/getpt.F [deleted file]
ISAJET/code/gettot.F [deleted file]
ISAJET/code/heavyx.F [deleted file]
ISAJET/code/hevolv.F [deleted file]
ISAJET/code/higgs.F [deleted file]
ISAJET/code/idanti.F [deleted file]
ISAJET/code/idgen.F [deleted file]
ISAJET/code/iframs.F [deleted file]
ISAJET/code/inisap.F [deleted file]
ISAJET/code/ipartns.F [deleted file]
ISAJET/code/ipjset.F [deleted file]
ISAJET/code/iprtns.F [deleted file]
ISAJET/code/irmov0.F [deleted file]
ISAJET/code/isabeg.F [deleted file]
ISAJET/code/isabg2.F [deleted file]
ISAJET/code/isaend.F [deleted file]
ISAJET/code/isaevt.F [deleted file]
ISAJET/code/isaini.F [deleted file]
ISAJET/code/isajet.F [deleted file]
ISAJET/code/isasrt.F [deleted file]
ISAJET/code/ispjet.F [deleted file]
ISAJET/code/istrad.F [deleted file]
ISAJET/code/iswdky.F [deleted file]
ISAJET/code/jetgen.F [deleted file]
ISAJET/code/kkgf1.F [deleted file]
ISAJET/code/kkgf2.F [deleted file]
ISAJET/code/kkgf3.F [deleted file]
ISAJET/code/label.F [deleted file]
ISAJET/code/lboost.F [deleted file]
ISAJET/code/logerr.F [deleted file]
ISAJET/code/logic.F [deleted file]
ISAJET/code/logmgm.F [deleted file]
ISAJET/code/logmgy.F [deleted file]
ISAJET/code/logmij.F [deleted file]
ISAJET/code/logp.F [deleted file]
ISAJET/code/logphi.F [deleted file]
ISAJET/code/logphw.F [deleted file]
ISAJET/code/logpt.F [deleted file]
ISAJET/code/logqm.F [deleted file]
ISAJET/code/logqt.F [deleted file]
ISAJET/code/logthw.F [deleted file]
ISAJET/code/logx.F [deleted file]
ISAJET/code/logxw.F [deleted file]
ISAJET/code/logyth.F [deleted file]
ISAJET/code/logyw.F [deleted file]
ISAJET/code/lstsq.F [deleted file]
ISAJET/code/mbias.F [deleted file]
ISAJET/code/mbset.F [deleted file]
ISAJET/code/mginit.F [deleted file]
ISAJET/code/muljet.F [deleted file]
ISAJET/code/nogood.F [deleted file]
ISAJET/code/ordecr.F [deleted file]
ISAJET/code/order.F [deleted file]
ISAJET/code/prtevt.F [deleted file]
ISAJET/code/prtlim.F [deleted file]
ISAJET/code/ptfun.F [deleted file]
ISAJET/code/qcdini.F [deleted file]
ISAJET/code/qcdint.F [deleted file]
ISAJET/code/qcdinz.F [deleted file]
ISAJET/code/qcdjet.F [deleted file]
ISAJET/code/qcdt.F [deleted file]
ISAJET/code/qcdz.F [deleted file]
ISAJET/code/qfunc.F [deleted file]
ISAJET/code/ranfgt.F [deleted file]
ISAJET/code/ranfmt.F [deleted file]
ISAJET/code/ranfst.F [deleted file]
ISAJET/code/readin.F [deleted file]
ISAJET/code/rejfrg.F [deleted file]
ISAJET/code/rejjet.F [deleted file]
ISAJET/code/rescal.F [deleted file]
ISAJET/code/reset.F [deleted file]
ISAJET/code/setcon.F [deleted file]
ISAJET/code/setdky.F [deleted file]
ISAJET/code/seth.F [deleted file]
ISAJET/code/sethss.F [deleted file]
ISAJET/code/setkkg.F [deleted file]
ISAJET/code/setnxt.F [deleted file]
ISAJET/code/settyp.F [deleted file]
ISAJET/code/setw.F [deleted file]
ISAJET/code/sigdy.F [deleted file]
ISAJET/code/sigdy2.F [deleted file]
ISAJET/code/sigee.F [deleted file]
ISAJET/code/sigfil.F [deleted file]
ISAJET/code/siggam.F [deleted file]
ISAJET/code/sigh.F [deleted file]
ISAJET/code/sigh2.F [deleted file]
ISAJET/code/sigh3.F [deleted file]
ISAJET/code/sighss.F [deleted file]
ISAJET/code/sigint.F [deleted file]
ISAJET/code/sigkkg.F [deleted file]
ISAJET/code/sigqcd.F [deleted file]
ISAJET/code/sigsse.F [deleted file]
ISAJET/code/sigssl.F [deleted file]
ISAJET/code/sigssy.F [deleted file]
ISAJET/code/sigssz.F [deleted file]
ISAJET/code/sigtc.F [deleted file]
ISAJET/code/sigtc2.F [deleted file]
ISAJET/code/sigtc3.F [deleted file]
ISAJET/code/sigwh.F [deleted file]
ISAJET/code/sigwhs.F [deleted file]
ISAJET/code/sigww.F [deleted file]
ISAJET/code/sigww2.F [deleted file]
ISAJET/code/smszg.F [deleted file]
ISAJET/code/spline.F [deleted file]
ISAJET/code/ssfel.F [deleted file]
ISAJET/code/ssgst.F [deleted file]
ISAJET/code/ssgt.F [deleted file]
ISAJET/code/struc.F [deleted file]
ISAJET/code/strucw.F [deleted file]
ISAJET/code/szjj1.F [deleted file]
ISAJET/code/szjj2.F [deleted file]
ISAJET/code/szjj3.F [deleted file]
ISAJET/code/szjj4.F [deleted file]
ISAJET/code/szjj5.F [deleted file]
ISAJET/code/szjj6.F [deleted file]
ISAJET/code/szjj7.F [deleted file]
ISAJET/code/timer.F [deleted file]
ISAJET/code/twojet.F [deleted file]
ISAJET/code/twokin.F [deleted file]
ISAJET/code/visaje.F [deleted file]
ISAJET/code/whiggs.F [deleted file]
ISAJET/code/wpair.F [deleted file]
ISAJET/code/wwkin.F [deleted file]
ISAJET/code/wwss.F [deleted file]
ISAJET/code/wwst.F [deleted file]
ISAJET/code/wwtt.F [deleted file]
ISAJET/code/wzss.F [deleted file]
ISAJET/code/wzst.F [deleted file]
ISAJET/code/wzsu.F [deleted file]
ISAJET/code/wztu.F [deleted file]
ISAJET/code/xwwww.F [deleted file]
ISAJET/code/xwwzz.F [deleted file]
ISAJET/code/xzzww.F [deleted file]
ISAJET/code/xzzzz.F [deleted file]
ISAJET/code/ygenj.F [deleted file]
ISAJET/code/zjj.F [deleted file]
ISAJET/code/zjj0.F [deleted file]
ISAJET/code/zjj1.F [deleted file]
ISAJET/code/zjj2.F [deleted file]
ISAJET/code/zjj3.F [deleted file]
ISAJET/code/zjj4.F [deleted file]
ISAJET/code/zjj5.F [deleted file]
ISAJET/code/zjj6.F [deleted file]
ISAJET/code/zjj7.F [deleted file]
ISAJET/code/zzall.F [deleted file]
ISAJET/code/zzstar.F [deleted file]
ISAJET/data/decay.cpp [deleted file]
ISAJET/doc/changes.doc [deleted file]
ISAJET/doc/decay.doc [deleted file]
ISAJET/doc/higher.doc [deleted file]
ISAJET/doc/ident.doc [deleted file]
ISAJET/doc/input.doc [deleted file]
ISAJET/doc/intro.doc [deleted file]
ISAJET/doc/isassdoc.doc [deleted file]
ISAJET/doc/main.doc [deleted file]
ISAJET/doc/output.doc [deleted file]
ISAJET/doc/patchy.doc [deleted file]
ISAJET/doc/physics.doc [deleted file]
ISAJET/doc/sample.doc [deleted file]
ISAJET/doc/susy.doc [deleted file]
ISAJET/doc/tape.doc [deleted file]
ISAJET/doc/ztext.doc [deleted file]
ISAJET/isadata/aldata.F [deleted file]
ISAJET/isajet/brembm.inc [deleted file]
ISAJET/isajet/calor.inc [deleted file]
ISAJET/isajet/const.inc [deleted file]
ISAJET/isajet/dkyss3.inc [deleted file]
ISAJET/isajet/dkytab.inc [deleted file]
ISAJET/isajet/dylim.inc [deleted file]
ISAJET/isajet/dypar.inc [deleted file]
ISAJET/isajet/eepar.inc [deleted file]
ISAJET/isajet/final.inc [deleted file]
ISAJET/isajet/force.inc [deleted file]
ISAJET/isajet/frame.inc [deleted file]
ISAJET/isajet/frgpar.inc [deleted file]
ISAJET/isajet/getjet.inc [deleted file]
ISAJET/isajet/hcon.inc [deleted file]
ISAJET/isajet/hcon1.inc [deleted file]
ISAJET/isajet/hcon2.inc [deleted file]
ISAJET/isajet/hepevt.inc [deleted file]
ISAJET/isajet/idrun.inc [deleted file]
ISAJET/isajet/isabnk.inc [deleted file]
ISAJET/isajet/isalnk.inc [deleted file]
ISAJET/isajet/isapw.inc [deleted file]
ISAJET/isajet/isaunt.inc [deleted file]
ISAJET/isajet/isloop.inc [deleted file]
ISAJET/isajet/ita.inc [deleted file]
ISAJET/isajet/itapes.inc [deleted file]
ISAJET/isajet/izisab.inc [deleted file]
ISAJET/isajet/izisac.inc [deleted file]
ISAJET/isajet/izisae.inc [deleted file]
ISAJET/isajet/izisaf.inc [deleted file]
ISAJET/isajet/izisaj.inc [deleted file]
ISAJET/isajet/izisal.inc [deleted file]
ISAJET/isajet/izisam.inc [deleted file]
ISAJET/isajet/izisaq.inc [deleted file]
ISAJET/isajet/iziscl.inc [deleted file]
ISAJET/isajet/iziscm.inc [deleted file]
ISAJET/isajet/izisjt.inc [deleted file]
ISAJET/isajet/izismr.inc [deleted file]
ISAJET/isajet/izisp1.inc [deleted file]
ISAJET/isajet/izisp2.inc [deleted file]
ISAJET/isajet/izisp3.inc [deleted file]
ISAJET/isajet/izisrc.inc [deleted file]
ISAJET/isajet/izisv1.inc [deleted file]
ISAJET/isajet/izisv2.inc [deleted file]
ISAJET/isajet/izpjet.inc [deleted file]
ISAJET/isajet/izpjhd.inc [deleted file]
ISAJET/isajet/izpjpt.inc [deleted file]
ISAJET/isajet/jetlim.inc [deleted file]
ISAJET/isajet/jetpar.inc [deleted file]
ISAJET/isajet/jetset.inc [deleted file]
ISAJET/isajet/jetsig.inc [deleted file]
ISAJET/isajet/jwork.inc [deleted file]
ISAJET/isajet/jwork2.inc [deleted file]
ISAJET/isajet/keys.inc [deleted file]
ISAJET/isajet/kkgrav.inc [deleted file]
ISAJET/isajet/l2cal.inc [deleted file]
ISAJET/isajet/l2dky.inc [deleted file]
ISAJET/isajet/l2getj.inc [deleted file]
ISAJET/isajet/l2jset.inc [deleted file]
ISAJET/isajet/l2part.inc [deleted file]
ISAJET/isajet/l2sigs.inc [deleted file]
ISAJET/isajet/l2zevl.inc [deleted file]
ISAJET/isajet/l2zout.inc [deleted file]
ISAJET/isajet/limevl.inc [deleted file]
ISAJET/isajet/listss.inc [deleted file]
ISAJET/isajet/lkpjet.inc [deleted file]
ISAJET/isajet/lstprt.inc [deleted file]
ISAJET/isajet/mbgen.inc [deleted file]
ISAJET/isajet/mbpar.inc [deleted file]
ISAJET/isajet/mgcoms.inc [deleted file]
ISAJET/isajet/mgkin.inc [deleted file]
ISAJET/isajet/mglims.inc [deleted file]
ISAJET/isajet/mgsigs.inc [deleted file]
ISAJET/isajet/myhist.inc [deleted file]
ISAJET/isajet/nodcay.inc [deleted file]
ISAJET/isajet/partcl.inc [deleted file]
ISAJET/isajet/pi.inc [deleted file]
ISAJET/isajet/pilot.h [deleted file]
ISAJET/isajet/pinits.inc [deleted file]
ISAJET/isajet/pjets.inc [deleted file]
ISAJET/isajet/primar.inc [deleted file]
ISAJET/isajet/prtout.inc [deleted file]
ISAJET/isajet/ptpar.inc [deleted file]
ISAJET/isajet/q1q2.inc [deleted file]
ISAJET/isajet/qcdpar.inc [deleted file]
ISAJET/isajet/qlmass.inc [deleted file]
ISAJET/isajet/qsave.inc [deleted file]
ISAJET/isajet/quest.inc [deleted file]
ISAJET/isajet/rectp.inc [deleted file]
ISAJET/isajet/seed.inc [deleted file]
ISAJET/isajet/ssinf.inc [deleted file]
ISAJET/isajet/sslun.inc [deleted file]
ISAJET/isajet/ssmode.inc [deleted file]
ISAJET/isajet/sspar.inc [deleted file]
ISAJET/isajet/sspols.inc [deleted file]
ISAJET/isajet/sssm.inc [deleted file]
ISAJET/isajet/sstmp.inc [deleted file]
ISAJET/isajet/sstype.inc [deleted file]
ISAJET/isajet/sugmg.inc [deleted file]
ISAJET/isajet/sugnu.inc [deleted file]
ISAJET/isajet/sugpas.inc [deleted file]
ISAJET/isajet/sugxin.inc [deleted file]
ISAJET/isajet/tcpar.inc [deleted file]
ISAJET/isajet/times.inc [deleted file]
ISAJET/isajet/totals.inc [deleted file]
ISAJET/isajet/types.inc [deleted file]
ISAJET/isajet/w50510.inc [deleted file]
ISAJET/isajet/w50517.inc [deleted file]
ISAJET/isajet/wcon.inc [deleted file]
ISAJET/isajet/wcon1.inc [deleted file]
ISAJET/isajet/wcon2.inc [deleted file]
ISAJET/isajet/wgen.inc [deleted file]
ISAJET/isajet/wsig.inc [deleted file]
ISAJET/isajet/wwpar.inc [deleted file]
ISAJET/isajet/wwpar1.inc [deleted file]
ISAJET/isajet/wwpar2.inc [deleted file]
ISAJET/isajet/wwsig.inc [deleted file]
ISAJET/isajet/xmssm.inc [deleted file]
ISAJET/isajet/zebcom.inc [deleted file]
ISAJET/isajet/zevel.inc [deleted file]
ISAJET/isajet/zlinka.inc [deleted file]
ISAJET/isajet/zvout.inc [deleted file]
ISAJET/isajetLinkDef.h [deleted file]
ISAJET/isarun/dialog.F [deleted file]
ISAJET/isarun/isaset.F [deleted file]
ISAJET/isasusy/ssalfs.F [deleted file]
ISAJET/isasusy/ssb0.F [deleted file]
ISAJET/isasusy/ssb1.F [deleted file]
ISAJET/isasusy/ssdhll.F [deleted file]
ISAJET/isasusy/ssdint.F [deleted file]
ISAJET/isasusy/ssdlam.F [deleted file]
ISAJET/isasusy/ssf0.F [deleted file]
ISAJET/isasusy/ssglbf.F [deleted file]
ISAJET/isasusy/ssgwq1.F [deleted file]
ISAJET/isasusy/ssgwq2.F [deleted file]
ISAJET/isasusy/ssgwt1.F [deleted file]
ISAJET/isasusy/ssgwt2.F [deleted file]
ISAJET/isasusy/ssgwt3.F [deleted file]
ISAJET/isasusy/ssgwt4.F [deleted file]
ISAJET/isasusy/ssgwt5.F [deleted file]
ISAJET/isasusy/ssgwt6.F [deleted file]
ISAJET/isasusy/ssgwt7.F [deleted file]
ISAJET/isasusy/ssgwt8.F [deleted file]
ISAJET/isasusy/ssgx1.F [deleted file]
ISAJET/isasusy/ssgx10.F [deleted file]
ISAJET/isasusy/ssgx11.F [deleted file]
ISAJET/isasusy/ssgx2.F [deleted file]
ISAJET/isasusy/ssgx3.F [deleted file]
ISAJET/isasusy/ssgx4.F [deleted file]
ISAJET/isasusy/ssgx5.F [deleted file]
ISAJET/isasusy/ssgx6.F [deleted file]
ISAJET/isasusy/ssgx7.F [deleted file]
ISAJET/isasusy/ssgx8.F [deleted file]
ISAJET/isasusy/ssgx9.F [deleted file]
ISAJET/isasusy/ssgzg1.F [deleted file]
ISAJET/isasusy/ssgzg2.F [deleted file]
ISAJET/isasusy/ssgzg3.F [deleted file]
ISAJET/isasusy/ssgzt.F [deleted file]
ISAJET/isasusy/sshcc.F [deleted file]
ISAJET/isasusy/sshff.F [deleted file]
ISAJET/isasusy/sshff1.F [deleted file]
ISAJET/isasusy/sshgl.F [deleted file]
ISAJET/isasusy/sshgm.F [deleted file]
ISAJET/isasusy/sshgm1.F [deleted file]
ISAJET/isasusy/sshhx.F [deleted file]
ISAJET/isasusy/sshibf.F [deleted file]
ISAJET/isasusy/sshnn.F [deleted file]
ISAJET/isasusy/sshsf.F [deleted file]
ISAJET/isasusy/sshww.F [deleted file]
ISAJET/isasusy/sshww1.F [deleted file]
ISAJET/isasusy/sshww2.F [deleted file]
ISAJET/isasusy/ssl1st.F [deleted file]
ISAJET/isasusy/sslpbf.F [deleted file]
ISAJET/isasusy/sslrt1.F [deleted file]
ISAJET/isasusy/ssmass.F [deleted file]
ISAJET/isasusy/ssme3.F [deleted file]
ISAJET/isasusy/ssmhc.F [deleted file]
ISAJET/isasusy/ssmhn.F [deleted file]
ISAJET/isasusy/ssmqcd.F [deleted file]
ISAJET/isasusy/ssmssm.F [deleted file]
ISAJET/isasusy/ssn1st.F [deleted file]
ISAJET/isasusy/ssnorm.F [deleted file]
ISAJET/isasusy/sspole.F [deleted file]
ISAJET/isasusy/ssqkbf.F [deleted file]
ISAJET/isasusy/sssave.F [deleted file]
ISAJET/isasusy/sssnws.F [deleted file]
ISAJET/isasusy/ssstbf.F [deleted file]
ISAJET/isasusy/sstest.F [deleted file]
ISAJET/isasusy/sstpbf.F [deleted file]
ISAJET/isasusy/sswibf.F [deleted file]
ISAJET/isasusy/sswwf1.F [deleted file]
ISAJET/isasusy/sswzbf.F [deleted file]
ISAJET/isasusy/sswzf1.F [deleted file]
ISAJET/isasusy/sswzf2.F [deleted file]
ISAJET/isasusy/sswzf3.F [deleted file]
ISAJET/isasusy/sswzf4.F [deleted file]
ISAJET/isasusy/sswzf5.F [deleted file]
ISAJET/isasusy/sswzf6.F [deleted file]
ISAJET/isasusy/sswzf7.F [deleted file]
ISAJET/isasusy/ssxint.F [deleted file]
ISAJET/isasusy/ssxlam.F [deleted file]
ISAJET/isasusy/sszhx.F [deleted file]
ISAJET/isasusy/sszibf.F [deleted file]
ISAJET/isasusy/sszwf1.F [deleted file]
ISAJET/isasusy/sszzf1.F [deleted file]
ISAJET/isasusy/sszzf2.F [deleted file]
ISAJET/isasusy/sszzf3.F [deleted file]
ISAJET/isasusy/sszzf4.F [deleted file]
ISAJET/isasusy/sszzf5.F [deleted file]
ISAJET/isasusy/sualfe.F [deleted file]
ISAJET/isasusy/sualfs.F [deleted file]
ISAJET/isasusy/sugeff.F [deleted file]
ISAJET/isasusy/sugfrz.F [deleted file]
ISAJET/isasusy/sugmas.F [deleted file]
ISAJET/isasusy/sugra.F [deleted file]
ISAJET/isasusy/sugrge.F [deleted file]
ISAJET/isasusy/surg06.F [deleted file]
ISAJET/isasusy/surg26.F [deleted file]
ISAJET/isatape/bufin.F [deleted file]
ISAJET/isatape/bufout.F [deleted file]
ISAJET/isatape/edit.F [deleted file]
ISAJET/isatape/isahep.F [deleted file]
ISAJET/isatape/isawbg.F [deleted file]
ISAJET/isatape/isawev.F [deleted file]
ISAJET/isatape/isawnd.F [deleted file]
ISAJET/isatape/itrans.F [deleted file]
ISAJET/isatape/movlev.F [deleted file]
ISAJET/isatape/prtlst.F [deleted file]
ISAJET/isatape/rdbeg.F [deleted file]
ISAJET/isatape/rdtape.F [deleted file]
ISAJET/isatape/rend.F [deleted file]
ISAJET/isatape/rgens.F [deleted file]
ISAJET/isatape/wgens.F [deleted file]
ISAJET/isatape/zerol.F [deleted file]
ISAJET/libisajet.pkg [deleted file]
ISAJET/main.c [deleted file]
ISAJET/openfile/openfile.F [deleted file]
ISAJET/pdfinit/pdfinit.F [deleted file]
ISAJET/utils/cern_lib/ddilog.F [deleted file]
ISAJET/utils/cern_lib/eisrs1.F [deleted file]
ISAJET/utils/cern_lib/rkstp.F [deleted file]
ISAJET/utils/cern_lib/sorttf.F [deleted file]
ISAJET/utils/cern_lib/tql2.F [deleted file]
ISAJET/utils/cern_lib/tred2.F [deleted file]
TIsajet/AliGenIsajet.cxx [deleted file]
TIsajet/AliGenIsajet.h [deleted file]
TIsajet/CMakeLists.txt [deleted file]
TIsajet/CMake_libTIsajet.txt [deleted file]
TIsajet/Icommon.h [deleted file]
TIsajet/TIsajet.cxx [deleted file]
TIsajet/TIsajet.h [deleted file]
TIsajet/TIsajetLinkDef.h [deleted file]
TIsajet/libTIsajet.pkg [deleted file]

diff --git a/ISAJET/AliIsajetRndm.cxx b/ISAJET/AliIsajetRndm.cxx
deleted file mode 100644 (file)
index 275abb9..0000000
+++ /dev/null
@@ -1,72 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-/* $Id$ */
-
-//-----------------------------------------------------------------------------
-//   Class: AliHerwigRndm
-//   Responsibilities: Interface to Root random number generator 
-//                     from Fortran (re-implements FINCTION RLU_HERWIG 
-//                     from HERWIG)
-//   Note: Since AliGenHerwig belongs to another module (THerwig) one cannot
-//         pass the ponter to the generator via static variable
-//   Collaborators: AliGenHerwig class
-//   Example:
-//
-//   root> AliGenHerwig *gener = new AliGenHerwig(-1);
-//   root> AliHerwigRndm::SetHerwigRandom(new TRandom3());
-//   root> AliHerwigRndm::GetHerwigRandom()->SetSeed(0);
-//   root> cout<<"Seed "<< AliHerwigRndm::GetHerwigRandom()->GetSeed() <<endl;
-//-----------------------------------------------------------------------------
-
-#include <TRandom.h>
-
-#include "AliIsajetRndm.h"
-
-TRandom * AliIsajetRndm::fgIsajetRandom=0;
-
-ClassImp(AliIsajetRndm)
-
-//_______________________________________________________________________
-void AliIsajetRndm::SetIsajetRandom(TRandom *ran) {
-  //
-  // Sets the pointer to an existing random numbers generator
-  //
-  if(ran) fgIsajetRandom=ran;
-  else fgIsajetRandom=gRandom;
-}
-
-//_______________________________________________________________________
-TRandom * AliIsajetRndm::GetIsajetRandom() {
-  //
-  // Retrieves the pointer to the random numbers generator
-  //
-  return fgIsajetRandom;
-}
-
-//_______________________________________________________________________
-# define ranf    ranf_
-
-extern "C" {
-  Double_t ranf() 
-  {
-    // Wrapper to FUNCTION HWR from HERWIG
-    // Uses static method to retrieve the pointer to the (C++) generator
-      Double_t r;
-      do r=AliIsajetRndm::GetIsajetRandom()->Rndm(); 
-      while(0 >= r || r >= 1);
-      return r;
-  }
-}
diff --git a/ISAJET/AliIsajetRndm.h b/ISAJET/AliIsajetRndm.h
deleted file mode 100644 (file)
index 5938797..0000000
+++ /dev/null
@@ -1,44 +0,0 @@
-#ifndef ALIISAJETRNDM_H
-#define ALIISAJETRNDM_H
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-/* $Id$ */
-
-#include <Rtypes.h>
-#include <TError.h>
-
-class TRandom;
-
-class AliIsajetRndm {
- public:
-  AliIsajetRndm() {
-    // Default constructor. The static data member is initialized
-    // in the implementation file
-  }
-  AliIsajetRndm(const AliIsajetRndm &/*rn*/) {
-    // Copy constructor: no copy allowed for the object
-    ::Fatal("Copy constructor","Not allowed\n");
-  }
-  virtual ~AliIsajetRndm() {
-    // Destructor
-    fgIsajetRandom=0;
-  }
-  AliIsajetRndm & operator=(const AliIsajetRndm &/*rn*/) {
-    // Assignment operator: no assignment allowed
-    ::Fatal("Assignment operator","Not allowed\n");
-    return (*this);
-  }
-  
-  static void SetIsajetRandom(TRandom *ran=0);
-  static TRandom * GetIsajetRandom();
-
-private:
-
-  static TRandom * fgIsajetRandom; //! pointer to the random number generator
-
-  ClassDef(AliIsajetRndm,0)  //Random Number generator wrapper (non persistent)
-};
-
-#endif 
-
diff --git a/ISAJET/CMakeLists.txt b/ISAJET/CMakeLists.txt
deleted file mode 100644 (file)
index 04224ad..0000000
+++ /dev/null
@@ -1,26 +0,0 @@
-# -*- mode: cmake -*-
-# Create a library called "lib<name>" which includes the source files given in
-# the array .
-# The extension is already found.  Any number of sources could be listed here.
-
-set(INCLUDE_DIRECTORIES
-${CMAKE_SOURCE_DIR}/ISAJET
-${CMAKE_SOURCE_DIR}/ISAJET/isajet
-${ROOT_INCLUDE_DIR}
-)
-
-include_directories( ${INCLUDE_DIRECTORIES})
-
-set(LINK_DIRECTORIES
-${ROOT_LIBRARY_DIR}
-) 
-
-link_directories( ${LINK_DIRECTORIES})
-
-SetModule()
-
-include(CMake_libisajet.txt)
-
-
-
-
diff --git a/ISAJET/CMake_libisajet.txt b/ISAJET/CMake_libisajet.txt
deleted file mode 100644 (file)
index 855ef09..0000000
+++ /dev/null
@@ -1,363 +0,0 @@
-# -*- mode: cmake -*-
-
-# Create a library called "lib<name>" which includes the source files given in
-# the array .
-# The extension is already found.  Any number of sources could be listed here.
-
-set(isajet_SRCS
-AliIsajetRndm.cxx
-)
-
-set(isajet_CSRCS
-main.c
-)
-
-set(isajet_FSRCS
-code/alqcd.F 
-code/amass.F 
-code/amgmw.F 
-code/charge.F 
-code/ctxc2i.F 
-code/ctxi2c.F 
-code/ctxin.F 
-code/ctxout.F 
-code/datime.F 
-code/dblpcm.F 
-code/dblvec.F 
-code/dboost.F 
-code/decay.F 
-code/decjet.F 
-code/decps1.F 
-code/decps2.F 
-code/decss3.F 
-code/dectau.F 
-code/decva.F 
-code/dhelas.F 
-code/dincgm.F 
-code/domssm.F 
-code/drllyn.F 
-code/ebeam.F 
-code/eebeg.F 
-code/eemax.F 
-code/elctrn.F 
-code/epf.F 
-code/estruc.F 
-code/evol01.F 
-code/evol02.F 
-code/evol03.F 
-code/evol05.F 
-code/evol06.F 
-code/evol07.F 
-code/evol11.F 
-code/evolms.F 
-code/evolve.F 
-code/fbrbm.F 
-code/flavor.F 
-code/fortop.F 
-code/frgjet.F 
-code/frgmnt.F 
-code/gamma.F 
-code/getpt.F 
-code/gettot.F 
-code/heavyx.F 
-code/hevolv.F 
-code/higgs.F 
-code/idanti.F 
-code/idgen.F 
-code/iframs.F 
-code/inisap.F 
-code/ipartns.F 
-code/ipjset.F 
-code/iprtns.F 
-code/irmov0.F 
-code/isabeg.F 
-code/isabg2.F 
-code/isaend.F 
-code/isaevt.F 
-code/isaini.F 
-code/isajet.F 
-code/isasrt.F 
-code/ispjet.F 
-code/istrad.F 
-code/iswdky.F 
-code/jetgen.F 
-code/kkgf1.F 
-code/kkgf2.F 
-code/kkgf3.F 
-code/label.F 
-code/lboost.F 
-code/logerr.F 
-code/logic.F 
-code/logmgm.F 
-code/logmgy.F 
-code/logmij.F 
-code/logp.F 
-code/logphi.F 
-code/logphw.F 
-code/logpt.F 
-code/logqm.F 
-code/logqt.F 
-code/logthw.F 
-code/logx.F 
-code/logxw.F 
-code/logyth.F 
-code/logyw.F 
-code/lstsq.F 
-code/mbias.F 
-code/mbset.F 
-code/mginit.F 
-code/muljet.F 
-code/nogood.F 
-code/ordecr.F 
-code/order.F 
-code/prtevt.F 
-code/prtlim.F 
-code/ptfun.F 
-code/qcdini.F 
-code/qcdint.F 
-code/qcdinz.F 
-code/qcdjet.F 
-code/qcdt.F 
-code/qcdz.F 
-code/qfunc.F 
-code/ranfgt.F 
-code/ranfmt.F 
-code/ranfst.F 
-code/readin.F 
-code/rejfrg.F 
-code/rejjet.F 
-code/rescal.F 
-code/reset.F 
-code/setcon.F 
-code/setdky.F 
-code/seth.F 
-code/sethss.F 
-code/setkkg.F 
-code/setnxt.F 
-code/settyp.F 
-code/setw.F 
-code/sigdy2.F 
-code/sigdy.F 
-code/sigee.F 
-code/sigfil.F 
-code/siggam.F 
-code/sigh2.F 
-code/sigh3.F 
-code/sigh.F 
-code/sighss.F 
-code/sigint.F 
-code/sigkkg.F 
-code/sigqcd.F 
-code/sigsse.F 
-code/sigssl.F 
-code/sigssy.F 
-code/sigssz.F 
-code/sigtc2.F 
-code/sigtc3.F 
-code/sigtc.F 
-code/sigwh.F 
-code/sigwhs.F 
-code/sigww2.F 
-code/sigww.F 
-code/smszg.F 
-code/spline.F 
-code/ssfel.F 
-code/ssgst.F 
-code/ssgt.F 
-code/struc.F 
-code/strucw.F 
-code/szjj1.F 
-code/szjj2.F 
-code/szjj3.F 
-code/szjj4.F 
-code/szjj5.F 
-code/szjj6.F 
-code/szjj7.F 
-code/timer.F 
-code/twojet.F 
-code/twokin.F 
-code/visaje.F 
-code/whiggs.F 
-code/wpair.F 
-code/wwkin.F 
-code/wwss.F 
-code/wwst.F 
-code/wwtt.F 
-code/wzss.F 
-code/wzst.F 
-code/wzsu.F 
-code/wztu.F 
-code/xwwww.F 
-code/xwwzz.F 
-code/xzzww.F 
-code/xzzzz.F 
-code/ygenj.F 
-code/zjj0.F 
-code/zjj1.F 
-code/zjj2.F 
-code/zjj3.F 
-code/zjj4.F 
-code/zjj5.F 
-code/zjj6.F 
-code/zjj7.F 
-code/zjj.F 
-code/zzall.F 
-code/zzstar.F 
-isasusy/ssalfs.F 
-isasusy/ssb0.F 
-isasusy/ssb1.F 
-isasusy/ssdhll.F 
-isasusy/ssdint.F 
-isasusy/ssdlam.F 
-isasusy/ssf0.F 
-isasusy/ssglbf.F 
-isasusy/ssgwq1.F 
-isasusy/ssgwq2.F 
-isasusy/ssgwt1.F 
-isasusy/ssgwt2.F 
-isasusy/ssgwt3.F 
-isasusy/ssgwt4.F 
-isasusy/ssgwt5.F 
-isasusy/ssgwt6.F 
-isasusy/ssgwt7.F 
-isasusy/ssgwt8.F 
-isasusy/ssgx10.F 
-isasusy/ssgx11.F 
-isasusy/ssgx1.F 
-isasusy/ssgx2.F 
-isasusy/ssgx3.F 
-isasusy/ssgx4.F 
-isasusy/ssgx5.F 
-isasusy/ssgx6.F 
-isasusy/ssgx7.F 
-isasusy/ssgx8.F 
-isasusy/ssgx9.F 
-isasusy/ssgzg1.F 
-isasusy/ssgzg2.F 
-isasusy/ssgzg3.F 
-isasusy/ssgzt.F 
-isasusy/sshcc.F 
-isasusy/sshff1.F 
-isasusy/sshff.F 
-isasusy/sshgl.F 
-isasusy/sshgm1.F 
-isasusy/sshgm.F 
-isasusy/sshhx.F 
-isasusy/sshibf.F 
-isasusy/sshnn.F 
-isasusy/sshsf.F 
-isasusy/sshww1.F 
-isasusy/sshww2.F 
-isasusy/sshww.F 
-isasusy/ssl1st.F 
-isasusy/sslpbf.F 
-isasusy/sslrt1.F 
-isasusy/ssmass.F 
-isasusy/ssme3.F 
-isasusy/ssmhc.F 
-isasusy/ssmhn.F 
-isasusy/ssmqcd.F 
-isasusy/ssmssm.F 
-isasusy/ssn1st.F 
-isasusy/ssnorm.F 
-isasusy/sspole.F 
-isasusy/ssqkbf.F 
-isasusy/sssave.F 
-isasusy/sssnws.F 
-isasusy/ssstbf.F 
-isasusy/sstest.F 
-isasusy/sstpbf.F 
-isasusy/sswibf.F 
-isasusy/sswwf1.F 
-isasusy/sswzbf.F 
-isasusy/sswzf1.F 
-isasusy/sswzf2.F 
-isasusy/sswzf3.F 
-isasusy/sswzf4.F 
-isasusy/sswzf5.F 
-isasusy/sswzf6.F 
-isasusy/sswzf7.F 
-isasusy/ssxint.F 
-isasusy/ssxlam.F 
-isasusy/sszhx.F 
-isasusy/sszibf.F 
-isasusy/sszwf1.F 
-isasusy/sszzf1.F 
-isasusy/sszzf2.F 
-isasusy/sszzf3.F 
-isasusy/sszzf4.F 
-isasusy/sszzf5.F 
-isasusy/sualfe.F 
-isasusy/sualfs.F 
-isasusy/sugeff.F 
-isasusy/sugfrz.F 
-isasusy/sugmas.F 
-isasusy/sugra.F 
-isasusy/sugrge.F 
-isasusy/surg06.F 
-isasusy/surg26.F 
-isatape/bufin.F 
-isatape/bufout.F 
-isatape/edit.F 
-isatape/isahep.F 
-isatape/isawbg.F 
-isatape/isawev.F 
-isatape/isawnd.F 
-isatape/itrans.F 
-isatape/movlev.F 
-isatape/prtlst.F 
-isatape/rdbeg.F 
-isatape/rdtape.F 
-isatape/rend.F 
-isatape/rgens.F 
-isatape/wgens.F 
-isatape/zerol.F 
-isadata/aldata.F 
-isarun/dialog.F 
-isarun/isaset.F 
-openfile/openfile.F 
-pdfinit/pdfinit.F 
-utils/cern_lib/ddilog.F 
-utils/cern_lib/eisrs1.F 
-utils/cern_lib/rkstp.F 
-utils/cern_lib/sorttf.F 
-utils/cern_lib/tql2.F 
-utils/cern_lib/tred2.F 
-)
-
-If(RULE_CHECKER_FOUND)
-  CHECK_RULES("${isajet_SRCS}" "${INCLUDE_DIRECTORIES}" isajet_RULES)
-endIf(RULE_CHECKER_FOUND)
-
-## fill list of header files from list of source files
-## by exchanging the file extension
-CHANGE_FILE_EXTENSION(*.cxx *.h isajet_HEADERS "${isajet_SRCS}")
-
-set(isajet_LINKDEF isajetLinkDef.h)
-set(isajet_DICTIONARY ${CMAKE_CURRENT_BINARY_DIR}/isajetDict.cxx) 
-#
-ROOT_GENERATE_DICTIONARY("${isajet_HEADERS}" "${isajet_LINKDEF}" "${isajet_DICTIONARY}" "${INCLUDE_DIRECTORIES}")
-#
-set(isajet_SRCS ${isajet_SRCS} ${isajet_DICTIONARY})
-
-add_library(isajet SHARED ${isajet_SRCS} ${isajet_CSRCS} ${isajet_FSRCS})
-target_link_libraries(isajet ${ROOT_LIBRARIES})
-set_target_properties(isajet PROPERTIES  ${ALIROOT_LIBRARY_PROPERTIES})
-SET_TARGET_PROPERTIES(isajet PROPERTIES LINKER_LANGUAGE CXX)
-
-################ install ###################
-install(TARGETS isajet DESTINATION ${ALIROOT_INSTALL_DIR}/lib)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/ISAJET/code/alqcd.F b/ISAJET/code/alqcd.F
deleted file mode 100644 (file)
index 08e4c29..0000000
+++ /dev/null
@@ -1,31 +0,0 @@
-#include "isajet/pilot.h"
-      FUNCTION ALQCD(Q2)  
-C-----------------------------------------------------------------------
-C     Strong coupling formula from page 201 of Barger and Phillips:
-C     (using ALQCD4 for 4 flavor Lambda)
-C-----------------------------------------------------------------------
-      REAL Q2,AS,TH5,TH6,PI,ALQCD4
-      LOGICAL FIRST
-      SAVE FIRST,PI,TH5,TH6,ALQCD4
-      DATA FIRST/.TRUE./
-C
-      IF(FIRST) THEN
-        PI=4.*ATAN(1.)
-        TH5=4*AMASS(5)**2
-        TH6=4*AMASS(6)**2
-        ALQCD4=0.177
-        FIRST=.FALSE.
-      ENDIF
-      IF (Q2.LE.TH5)THEN
-        AS=12*PI/(25*LOG(Q2/ALQCD4**2))
-      ELSE IF(Q2.GT.TH5.AND.Q2.LE.TH6) THEN
-        AS=25*LOG(Q2/ALQCD4**2)-2*LOG(Q2/TH5)
-        AS=12*PI/AS
-      ELSEIF(Q2.GT.TH6)THEN
-        AS=25*LOG(Q2/ALQCD4**2)
-        AS=AS-2*(LOG(Q2/TH5)+LOG(Q2/TH6))
-        AS=12*PI/AS
-      ENDIF
-      ALQCD=AS
-      RETURN
-      END
diff --git a/ISAJET/code/amass.F b/ISAJET/code/amass.F
deleted file mode 100644 (file)
index 1584440..0000000
+++ /dev/null
@@ -1,136 +0,0 @@
-#include "isajet/pilot.h"
-      FUNCTION AMASS(ID)
-C
-C          Returns the mass of the particle with IDENT code ID.
-C          Quark-based IDENT code.
-C          Ver 7.10: Update masses and split B baryon degeneracy.
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-#include "isajet/itapes.inc"
-#include "isajet/qlmass.inc"
-      INTEGER ID
-      REAL AMASS
-      REAL AMMES0(10),AMMES1(10),AMBAR0(30),AMBAR1(30)
-      INTEGER IFL1,IFL2,IFL3,JSPIN,INDEX,IFL1A,IFL2A,IFL3A,IDA
-C
-C          0- meson mass table
-C          pi0, pi+, eta, k+, k0, etap, ad0, d-, ds-, etac
-C
-      DATA AMMES0/.13496,.13957,.54745,.49367,.49767,.95775,1.8645
-     $,1.8693,1.9688,2.9788/
-C
-C          1- meson mass table
-C          rho0, rho+, omega, k*+, k*0, phi, ad*0, d*-, d*s-, jpsi
-C
-      DATA AMMES1/.7681,.7681,.78195,.89159,.89610,1.0194,2.0071
-     $,2.0101,2.1103,3.0969/
-C
-C          1/2+ baryon mass table
-C          x,p,n,-,-,s+,s0,s-,l,xi0,xi-,x,x,x
-C          sc++,sc+,sc0,lc+,usc.,dsc.,ssc.,sdc.,suc.,ucc.,dcc.,scc.
-C
-      DATA AMBAR0/-1.,.93828,.93957,2*-1.,1.1894,1.1925,1.1974
-     $,1.1156,1.3149,1.3213,3*-1.,2.4527,2.4529,2.4525,2.2849
-     $,2.50,2.50,2.60,2.40,2.40,3.55,3.55,3.70,4*-1./
-C
-C          3/2+ baryon mass table
-C          dl++,dl+,dl0,dl-,-,s*+,s*0,s*-,x,xi*0,xi*-,om-,x,x
-C          uuc*,udc*,ddc*,x,usc*,dsc*,ssc*,x,x,,ucc*,dcc*,scc*,ccc*
-C
-      DATA AMBAR1/1.232,1.232,1.232,1.232,-1.,1.3823,1.3820
-     $,1.3875,-1.,1.5318,1.5350,1.6722,2*-1.
-     $,2.63,2.63,2.63,-1.,2.70,2.70,2.80,2*-1.,3.75,3.75
-     $,3.90,4.80,3*-1./
-C
-C          Entry
-C
-      AMASS=-1.
-      CALL FLAVOR(ID,IFL1,IFL2,IFL3,JSPIN,INDEX)
-      IDA=IABS(ID)
-      IFL1A=IABS(IFL1)
-      IFL2A=IABS(IFL2)
-      IFL3A=IABS(IFL3)
-      IF(IDA.GT.10000.OR.JSPIN.GT.1) GO TO 500
-C
-C          Diquarks
-C
-      IF(ID.NE.0.AND.MOD(ID,100).EQ.0) THEN
-        AMASS=AMLEP(IFL1A)+AMLEP(IFL2A)
-C
-C          b and t particles. Only a few b masses are known, but we
-C          guess a few others to make sure decays are allowed:
-C
-      ELSEIF(IFL3A.GT.4) THEN
-        IF(IDA.EQ.150.OR.IDA.EQ.250) THEN
-          AMASS=5.2786
-        ELSEIF(IDA.EQ.151.OR.IDA.EQ.251) THEN
-          AMASS=5.3246
-        ELSEIF(IDA.EQ.350) THEN
-          AMASS=5.3693
-        ELSEIF(IDA.EQ.351) THEN
-          AMASS=5.3693+0.04
-        ELSEIF(IDA.EQ.2150) THEN
-          AMASS=5.641
-        ELSEIF(IDA.EQ.1150.OR.IDA.EQ.1250.OR.IDA.EQ.2250) THEN
-          AMASS=5.641+0.171
-        ELSEIF(IDA.EQ.2151) THEN
-          AMASS=5.641+.04
-        ELSEIF(IDA.EQ.1151.OR.IDA.EQ.1251.OR.IDA.EQ.2251) THEN
-          AMASS=5.641+0.171+0.04
-        ELSE
-          AMASS=AMLEP(IFL2A)+AMLEP(IFL3A)-.03+.04*JSPIN
-          IF(IFL1.NE.0) AMASS=AMASS+AMLEP(IFL1A)
-        ENDIF
-C
-C          Quarks and leptons
-C
-      ELSEIF(IFL2.EQ.0) THEN
-        AMASS=AMLEP(INDEX)
-C
-C          Mesons
-C
-      ELSEIF(IFL1.EQ.0) THEN
-        INDEX=INDEX-36*JSPIN-NQLEP
-        INDEX=INDEX-13
-        AMASS=(1-JSPIN)*AMMES0(INDEX)+JSPIN*AMMES1(INDEX)
-C
-C          Baryons
-C
-      ELSE
-        INDEX=INDEX-109*JSPIN-36*NMES-NQLEP
-        INDEX=INDEX-13
-        AMASS=(1-JSPIN)*AMBAR0(INDEX)+JSPIN*AMBAR1(INDEX)
-      ENDIF
-      RETURN
-C
-C          Special hadrons - used only in B decays
-C
-500   IF(IDA.EQ.10121.OR.IDA.EQ.10111) THEN
-        AMASS=1.230
-      ELSEIF(IDA.EQ.10131.OR.IDA.EQ.10231) THEN
-        AMASS=1.273
-      ELSEIF(IDA.EQ.30131.OR.IDA.EQ.30231) THEN
-        AMASS=1.412
-      ELSEIF(IDA.EQ.132) THEN
-        AMASS=1.4254
-      ELSEIF(IDA.EQ.232) THEN
-        AMASS=1.4324
-      ELSEIF(IDA.EQ.10110) THEN
-        AMASS=0.980+0.020
-      ELSEIF(IDA.EQ.112) THEN
-        AMASS=1.275
-      ELSEIF(IDA.EQ.10441) THEN
-        AMASS=3.686
-      ELSEIF(IDA.EQ.20440) THEN
-        AMASS=3.4151
-      ELSEIF(IDA.EQ.20441) THEN
-        AMASS=3.51053
-      ELSEIF(IDA.EQ.20442) THEN
-        AMASS=3.56617
-      ELSE
-        AMASS=0
-      ENDIF
-      RETURN
-      END
diff --git a/ISAJET/code/amgmw.F b/ISAJET/code/amgmw.F
deleted file mode 100644 (file)
index 479278f..0000000
+++ /dev/null
@@ -1,36 +0,0 @@
-#include "isajet/pilot.h"
-      FUNCTION AMGMW(I,J)
-C
-C          Get masses and widths from ISAJET commons for MadGraph
-C          I = particle IDENT
-C          J = 1 for mass
-C            = 2 for width
-C            = 3 for sin^2(theta)
-C          Needed to avoid common block name clashes with MadGraph
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-#include "isajet/itapes.inc"
-#include "isajet/wcon.inc"
-#include "isajet/hcon.inc"
-#include "isajet/sstype.inc"
-      INTEGER I,J
-      REAL AMGMW,AMASS
-C
-      IF(J.EQ.1) THEN
-        AMGMW=AMASS(I)
-      ELSEIF(J.EQ.2.AND.I.EQ.IDW) THEN
-        AMGMW=WGAM(2)
-      ELSEIF(J.EQ.2.AND.I.EQ.IDZ) THEN
-        AMGMW=WGAM(4)
-      ELSEIF(J.EQ.2.AND.I.EQ.IDH) THEN
-        AMGMW=HGAM
-      ELSEIF(J.EQ.3.AND.I.EQ.1) THEN
-        AMGMW=SIN2W
-      ELSE
-        WRITE(ITLIS,*) 'ERROR IN AMGMW: I,J =',I,J
-        STOP99
-      ENDIF
-      RETURN
-      END
diff --git a/ISAJET/code/charge.F b/ISAJET/code/charge.F
deleted file mode 100644 (file)
index 951333c..0000000
+++ /dev/null
@@ -1,30 +0,0 @@
-#include "isajet/pilot.h"
-      FUNCTION CHARGE(ID)
-C
-C          COMPUTE CHARGE OF PARTICLE WITH IDENT CODE ID
-C          ICHRG MUST BE DIMENSIONED NQLEP+13
-C
-#include "isajet/itapes.inc"
-      DIMENSION ICHRG(75),IFL(3)
-C          3 * charge
-      DATA ICHRG/0
-     $,2,-1,-1,2,-1,2,-1,2,0,0, 0,-3,0,-3,0,-3,0,-3,0,0,0
-     $,2,-1,-1,2,-1,2,-1,2,0,0, 0,-3,0,-3,0,-3,0,-3,3,0
-     $,2,-1,-1,2,-1,2,-1,2,3,0, 0,-3,0,-3,0,-3,0,-3,3,0
-     $,3,0,0,0,0,0,3,3,6,6,0,0,0/
-C
-      IDABS=IABS(ID)
-      CALL FLAVOR(ID,IFL(1),IFL(2),IFL(3),JSPIN,INDEX)
-      IF(IDABS.LT.100) GO TO 200
-C
-      ISUM=0
-      DO 100 I=1,3
-        ISUM=ISUM+ICHRG(IABS(IFL(I))+1)*ISIGN(1,IFL(I))
-  100 CONTINUE
-      CHARGE=ISUM/3.
-      RETURN
-C
-200   CHARGE=ICHRG(INDEX+1)*ISIGN(1,ID)
-      CHARGE=CHARGE/3.
-      RETURN
-      END
diff --git a/ISAJET/code/ctxc2i.F b/ISAJET/code/ctxc2i.F
deleted file mode 100644 (file)
index add374e..0000000
+++ /dev/null
@@ -1,17 +0,0 @@
-#include "isajet/pilot.h"
-      SUBROUTINE CTXC2I(CVAL,IVAL,NSIZE)
-C-----------------------------------------------------------------------
-C          Convert character variable CVAL to integer array IVAL
-C-----------------------------------------------------------------------
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      CHARACTER*(*) CVAL
-      INTEGER I,NSIZE
-      INTEGER IVAL(NSIZE)
-C
-      DO 100 I=1,NSIZE
-100   IVAL(I)=ICHAR(CVAL(I:I))
-C
-      RETURN
-      END
diff --git a/ISAJET/code/ctxi2c.F b/ISAJET/code/ctxi2c.F
deleted file mode 100644 (file)
index 741a74d..0000000
+++ /dev/null
@@ -1,17 +0,0 @@
-#include "isajet/pilot.h"
-      SUBROUTINE CTXI2C(IVAL,CVAL,NSIZE)
-C-----------------------------------------------------------------------
-C          Convert integer array IVAL to character variable CVAL
-C-----------------------------------------------------------------------
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      CHARACTER*(*) CVAL
-      INTEGER I,NSIZE
-      INTEGER IVAL(NSIZE)
-C
-      DO 100 I=1,NSIZE
-100   CVAL(I:I)=CHAR(IVAL(I))
-C
-      RETURN
-      END
diff --git a/ISAJET/code/ctxin.F b/ISAJET/code/ctxin.F
deleted file mode 100644 (file)
index 47a199d..0000000
+++ /dev/null
@@ -1,201 +0,0 @@
-#include "isajet/pilot.h"
-      SUBROUTINE CTXIN(NVC,VC,MXVC)
-C-----------------------------------------------------------------------
-C  Purpose:
-C          Restore the context for an ISAJET job:
-C          Restore NVC words of VC all common blocks NOT associated only
-C          with a single event. Call CTXOUT and this to generate mixed
-C          events.
-C          PARAMETER (MXVC=20000)
-C          REAL    VC(MXVC)
-C          ...
-C          CALL CTXIN(NVC,VC,MXVC)
-C
-C          Note that the MSSM common blocks are not saved, so different
-C          SUSY runs cannot be mixed.
-C
-C          Ver. 7.02: Equivalenced dummy variables to avoid mixed 
-C                     arguments in MOVLEV or multiple EQUIVALENCEd
-C                     arguments to CTXIN/CTXOUT.
-C
-C  Author:
-C          F.E. Paige, April 1992     
-C-----------------------------------------------------------------------
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-#include "isajet/dkytab.inc"
-#include "isajet/dylim.inc"
-#include "isajet/dypar.inc"
-#include "isajet/eepar.inc"
-#include "isajet/final.inc"
-#include "isajet/force.inc"
-#include "isajet/frgpar.inc"
-#include "isajet/hcon.inc"
-#include "isajet/idrun.inc"
-#include "isajet/isloop.inc"
-#include "isajet/itapes.inc"
-#include "isajet/jetlim.inc"
-#include "isajet/keys.inc"
-#include "isajet/limevl.inc"
-#include "isajet/lstprt.inc"
-#include "isajet/mbgen.inc"
-#include "isajet/mbpar.inc"
-#include "isajet/nodcay.inc"
-#include "isajet/primar.inc"
-#include "isajet/prtout.inc"
-#include "isajet/ptpar.inc"
-#include "isajet/q1q2.inc"
-#include "isajet/qcdpar.inc"
-#include "isajet/qlmass.inc"
-#include "isajet/tcpar.inc"
-#include "isajet/times.inc"
-#include "isajet/totals.inc"
-#include "isajet/types.inc"
-#include "isajet/wcon.inc"
-C
-      INTEGER NVC,MXVC,NC,NN,I
-      REAL VC(MXVC)
-      CHARACTER*8 CLIST(290)
-      EQUIVALENCE (CLIST(1),PARTYP(1))
-C
-C          Dummy real variables for integers
-      REAL VLOOK(MXLOOK+6*MXDKY)
-      EQUIVALENCE (VLOOK(1),LOOK(1))
-      REAL VNKINF(5)
-      EQUIVALENCE (VNKINF(1),NKINF)
-      REAL VFORCE(9*MXFORC+1)
-      EQUIVALENCE (VFORCE(1),NFORCE)
-      REAL VIDVER(5)
-      EQUIVALENCE (VIDVER(1),IDVER)
-      REAL VEVOLV(4)
-      EQUIVALENCE (VEVOLV(1),NEVOLV)
-      REAL VITDKY(4)
-      EQUIVALENCE (VITDKY(1),ITDKY)
-      REAL VIKEYS(12)
-      EQUIVALENCE (VIKEYS(1),IKEYS)
-      REAL VSTPRT
-      EQUIVALENCE (VSTPRT,LSTPRT)
-      REAL VNJET(9)
-      EQUIVALENCE (VNJET(1),NJET)
-      REAL VEVPRT(2)
-      EQUIVALENCE (VEVPRT(1),NEVPRT)
-      REAL VKINPT(5)
-      EQUIVALENCE (VKINPT(1),NKINPT)
-      REAL VLOC(100)
-      EQUIVALENCE (VLOC(1),LOC(1))
-C          Dummy real variables for logicals
-      REAL VFLW(13)
-      EQUIVALENCE (VFLW(1),FLW)
-      REAL VNODCY(6)
-      EQUIVALENCE (VNODCY(1),NODCAY)
-      REAL VGOQ(3*MXGOQ+135)
-      EQUIVALENCE (VGOQ(1),GOQ(1,1))
-C
-      NC=0
-C          DKYTAB
-      NN=MXLOOK+6*MXDKY
-      CALL MOVLEV(VC(NC+1),VLOOK(1),NN)
-      NC=NC+NN
-C          DYLIM
-      CALL MOVLEV(VC(NC+1),QMIN,24)
-      NC=NC+24
-C          DYPAR
-      CALL MOVLEV(VC(NC+1),VFLW(1),13)
-      NC=NC+13
-C          EEPAR
-      CALL MOVLEV(VC(NC+1),SGMXEE,1)
-      NC=NC+1
-C          FINAL
-      CALL MOVLEV(VC(NC+1),VNKINF(1),5)
-      NC=NC+5
-C          FORCE
-      NN=9*MXFORC+1
-      CALL MOVLEV(VC(NC+1),VFORCE(1),NN)
-      NC=NC+NN
-C          FRGPAR
-      CALL MOVLEV(VC(NC+1),PUD,41)
-      NC=NC+41
-C          HCON
-      CALL MOVLEV(VC(NC+1),HMASS,69)
-      NC=NC+69
-C          IDRUN
-      CALL MOVLEV(VC(NC+1),VIDVER(1),5)
-      NC=NC+5
-C          ISLOOP
-      CALL MOVLEV(VC(NC+1),VEVOLV(1),4)
-      NC=NC+4
-C          ITAPES
-      CALL MOVLEV(VC(NC+1),VITDKY(1),4)
-      NC=NC+4
-C          JETLIM
-      CALL MOVLEV(VC(NC+1),PMIN(1),72)
-      NC=NC+72
-C          KEYS
-      CALL MOVLEV(VC(NC+1),VIKEYS(1),12)
-      NC=NC+12
-      CALL CTXI2C(VC(NC+1),REAC,8)
-      NC=NC+8
-C          LIMEVL
-      CALL MOVLEV(VC(NC+1),ETTHRS,3)
-      NC=NC+3
-C          LSTPRT
-      CALL MOVLEV(VC(NC+1),VSTPRT,1)
-      NC=NC+1
-C          MBGEN
-      NN=4*LIMPOM+8
-      CALL MOVLEV(VC(NC+1),POMWT(1),NN)
-      NC=NC+NN
-C          MBPAR
-      CALL MOVLEV(VC(NC+1),PUD0,19)
-      NC=NC+19
-C          NODCAY
-      CALL MOVLEV(VC(NC+1),VNODCY(1),6)
-      NC=NC+6
-C          PRIMAR
-      CALL MOVLEV(VC(NC+1),VNJET(1),9)
-      NC=NC+9
-C          PRTOUT
-      CALL MOVLEV(VC(NC+1),VEVPRT(1),2)
-      NC=NC+2
-C          PTPAR
-      CALL MOVLEV(VC(NC+1),PTFUN1,6)
-      NC=NC+6
-C          Q1Q2
-      CALL MOVLEV(VC(NC+1),VGOQ(1),3*MXGOQ+135)
-      NC=NC+3*MXGOQ+135
-C          QCDPAR
-      CALL MOVLEV(VC(NC+1),ALAM,4)
-      NC=NC+4
-C          QLMASS
-      CALL MOVLEV(VC(NC+1),AMLEP(1),55)
-      NC=NC+55
-C          TCPAR
-      CALL MOVLEV(VC(NC+1),TCMRHO,2)
-      NC=NC+2
-C          TIMES
-      CALL MOVLEV(VC(NC+1),TIME1,2)
-      NC=NC+2
-C          TOTALS
-      CALL MOVLEV(VC(NC+1),VKINPT(1),5)
-      NC=NC+5
-C          TYPES
-      CALL MOVLEV(VC(NC+1),VLOC(1),100)
-      NC=NC+100
-      DO 100 I=1,290
-        CALL CTXI2C(VC(NC+1),CLIST(I),8)
-        NC=NC+8
-100   CONTINUE
-C          WCON
-#if defined(CERNLIB_SINGLE)
-      NN=514
-#endif
-#if defined(CERNLIB_DOUBLE)
-      NN=514+97
-#endif
-      CALL MOVLEV(VC(NC+1),SIN2W,NN)
-      NC=NC+NN
-C
-      NVC=NC
-      RETURN
-      END
diff --git a/ISAJET/code/ctxout.F b/ISAJET/code/ctxout.F
deleted file mode 100644 (file)
index 4a2c0d7..0000000
+++ /dev/null
@@ -1,207 +0,0 @@
-#include "isajet/pilot.h"
-      SUBROUTINE CTXOUT(NVC,VC,MXVC)
-C-----------------------------------------------------------------------
-C  Purpose:
-C          Save the context for an ISAJET job:
-C          Save in NVC words of VC all common blocks NOT associated only
-C          with a single event. Call this and CTXIN to generate mixed
-C          events.
-C          PARAMETER (MXVC=20000)
-C          REAL    VC(MXVC)
-C          ...
-C          CALL CTXIN(NVC,VC,MXVC)
-C
-C          Note that the MSSM common blocks are not saved, so different
-C          SUSY runs cannot be mixed.
-C
-C          Ver. 7.02: Equivalenced dummy variables to avoid mixed 
-C                     arguments in MOVLEV or multiple EQUIVALENCEd
-C                     arguments to CTXIN/CTXOUT.
-C
-C  Author:
-C          F.E. Paige, April 1992     
-C-----------------------------------------------------------------------
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-#include "isajet/dkytab.inc"
-#include "isajet/dylim.inc"
-#include "isajet/dypar.inc"
-#include "isajet/eepar.inc"
-#include "isajet/final.inc"
-#include "isajet/force.inc"
-#include "isajet/frgpar.inc"
-#include "isajet/hcon.inc"
-#include "isajet/idrun.inc"
-#include "isajet/isloop.inc"
-#include "isajet/itapes.inc"
-#include "isajet/jetlim.inc"
-#include "isajet/keys.inc"
-#include "isajet/limevl.inc"
-#include "isajet/lstprt.inc"
-#include "isajet/mbgen.inc"
-#include "isajet/mbpar.inc"
-#include "isajet/nodcay.inc"
-#include "isajet/primar.inc"
-#include "isajet/prtout.inc"
-#include "isajet/ptpar.inc"
-#include "isajet/q1q2.inc"
-#include "isajet/qcdpar.inc"
-#include "isajet/qlmass.inc"
-#include "isajet/tcpar.inc"
-#include "isajet/times.inc"
-#include "isajet/totals.inc"
-#include "isajet/types.inc"
-#include "isajet/wcon.inc"
-C
-      INTEGER NVC,MXVC,NC,NN,I
-      REAL VC(MXVC)
-      CHARACTER*8 CLIST(290)
-      EQUIVALENCE (CLIST(1),PARTYP(1))
-C
-C          Dummy real variables for integers
-      REAL VLOOK(MXLOOK+6*MXDKY)
-      EQUIVALENCE (VLOOK(1),LOOK(1))
-      REAL VNKINF(5)
-      EQUIVALENCE (VNKINF(1),NKINF)
-      REAL VFORCE(9*MXFORC+1)
-      EQUIVALENCE (VFORCE(1),NFORCE)
-      REAL VIDVER(5)
-      EQUIVALENCE (VIDVER(1),IDVER)
-      REAL VEVOLV(4)
-      EQUIVALENCE (VEVOLV(1),NEVOLV)
-      REAL VITDKY(4)
-      EQUIVALENCE (VITDKY(1),ITDKY)
-      REAL VIKEYS(12)
-      EQUIVALENCE (VIKEYS(1),IKEYS)
-      REAL VSTPRT
-      EQUIVALENCE (VSTPRT,LSTPRT)
-      REAL VNJET(9)
-      EQUIVALENCE (VNJET(1),NJET)
-      REAL VEVPRT(2)
-      EQUIVALENCE (VEVPRT(1),NEVPRT)
-      REAL VKINPT(5)
-      EQUIVALENCE (VKINPT(1),NKINPT)
-      REAL VLOC(100)
-      EQUIVALENCE (VLOC(1),LOC(1))
-C          Dummy real variables for logicals
-      REAL VFLW(13)
-      EQUIVALENCE (VFLW(1),FLW)
-      REAL VNODCY(6)
-      EQUIVALENCE (VNODCY(1),NODCAY)
-      REAL VGOQ(3*MXGOQ+135)
-      EQUIVALENCE (VGOQ(1),GOQ(1,1))
-C
-      NC=0
-C          DKYTAB
-      NN=MXLOOK+6*MXDKY
-      CALL MOVLEV(VLOOK(1),VC(NC+1),NN)
-      NC=NC+NN
-C          DYLIM
-      CALL MOVLEV(QMIN,VC(NC+1),24)
-      NC=NC+24
-C          DYPAR
-      CALL MOVLEV(VFLW(1),VC(NC+1),13)
-      NC=NC+13
-C          EEPAR
-      CALL MOVLEV(SGMXEE,VC(NC+1),1)
-      NC=NC+1
-C          FINAL
-      CALL MOVLEV(VNKINF(1),VC(NC+1),5)
-      NC=NC+5
-C          FORCE
-      NN=9*MXFORC+1
-      CALL MOVLEV(VFORCE(1),VC(NC+1),NN)
-      NC=NC+NN
-C          FRGPAR
-      CALL MOVLEV(PUD,VC(NC+1),41)
-      NC=NC+41
-C          HCON
-      CALL MOVLEV(HMASS,VC(NC+1),69)
-      NC=NC+69
-C          IDRUN
-      CALL MOVLEV(VIDVER(1),VC(NC+1),5)
-      NC=NC+5
-C          ISLOOP
-      CALL MOVLEV(VEVOLV(1),VC(NC+1),4)
-      NC=NC+4
-C          ITAPES
-      CALL MOVLEV(VITDKY(1),VC(NC+1),4)
-      NC=NC+4
-C          JETLIM
-      CALL MOVLEV(PMIN(1),VC(NC+1),72)
-      NC=NC+72
-C          KEYS
-      CALL MOVLEV(VIKEYS(1),VC(NC+1),12)
-      NC=NC+12
-      CALL CTXC2I(REAC,VC(NC+1),8)
-      NC=NC+8
-C          LIMEVL
-      CALL MOVLEV(ETTHRS,VC(NC+1),3)
-      NC=NC+3
-C          LSTPRT
-      CALL MOVLEV(VSTPRT,VC(NC+1),1)
-      NC=NC+1
-C          MBGEN
-      NN=4*LIMPOM+8
-      CALL MOVLEV(POMWT(1),VC(NC+1),NN)
-      NC=NC+NN
-C          MBPAR
-      CALL MOVLEV(PUD0,VC(NC+1),19)
-      NC=NC+19
-C          NODCAY
-      CALL MOVLEV(VNODCY(1),VC(NC+1),6)
-      NC=NC+6
-C          PRIMAR
-      CALL MOVLEV(VNJET(1),VC(NC+1),9)
-      NC=NC+9
-C          PRTOUT
-      CALL MOVLEV(VEVPRT(1),VC(NC+1),2)
-      NC=NC+2
-C          PTPAR
-      CALL MOVLEV(PTFUN1,VC(NC+1),6)
-      NC=NC+6
-C          Q1Q2
-      CALL MOVLEV(VGOQ(1),VC(NC+1),3*MXGOQ+135)
-      NC=NC+3*MXGOQ+135
-C          QCDPAR
-      CALL MOVLEV(ALAM,VC(NC+1),4)
-      NC=NC+4
-C          QLMASS
-      CALL MOVLEV(AMLEP(1),VC(NC+1),55)
-      NC=NC+55
-C          TCPAR
-      CALL MOVLEV(TCMRHO,VC(NC+1),2)
-      NC=NC+2
-C          TIMES
-      CALL MOVLEV(TIME1,VC(NC+1),2)
-      NC=NC+2
-C          TOTALS
-      CALL MOVLEV(VKINPT(1),VC(NC+1),5)
-      NC=NC+5
-C          TYPES
-      CALL MOVLEV(VLOC(1),VC(NC+1),100)
-      NC=NC+100
-      DO 100 I=1,290
-        CALL CTXC2I(CLIST(I),VC(NC+1),8)
-        NC=NC+8
-100   CONTINUE
-C          WCON
-#if defined(CERNLIB_SINGLE)
-      NN=514
-#endif
-#if defined(CERNLIB_DOUBLE)
-      NN=514+97
-#endif
-      CALL MOVLEV(SIN2W,VC(NC+1),NN)
-      NC=NC+NN
-C
-      IF(NC.LE.MXVC) THEN
-        NVC=NC
-        RETURN
-      ELSE
-        WRITE(ITLIS,9000) NC
-9000    FORMAT(//' ERROR IN CTXOUT, NC = ',I5)
-        STOP99
-      ENDIF
-      END
diff --git a/ISAJET/code/datime.F b/ISAJET/code/datime.F
deleted file mode 100644 (file)
index 90395dd..0000000
+++ /dev/null
@@ -1,14 +0,0 @@
-#include "isajet/pilot.h"
-#if (defined(CERNLIB_VAX))&&(defined(CERNLIB_NOCERN))
-      SUBROUTINE DATIME(ID,IT)
-C          CALL VAX DATE AND TIME.
-#include "isajet/itapes.inc"
-      CHARACTER*8 BUF
-      CALL IDATE(IMON,IDAY,IYR)
-      CALL TIME(BUF)
-      ID=10000*IYR+100*IMON+IDAY
-      READ(BUF,'(I2,1X,I2,1X,I2)') K1,K2,K3
-      IT=10000*K1+100*K2+K3
-      RETURN
-      END
-#endif
diff --git a/ISAJET/code/dblpcm.F b/ISAJET/code/dblpcm.F
deleted file mode 100644 (file)
index b3f4a4d..0000000
+++ /dev/null
@@ -1,19 +0,0 @@
-#include "isajet/pilot.h"
-      FUNCTION DBLPCM(A,B,C)
-C          Calculate com momentum for A-->B+C with double precision.
-C          Needed to fix bug on 32-bit machines at high energy.
-C          Ver. 7.27: Rewrite order and then take abs value to be sure.
-#include "isajet/itapes.inc"
-#if defined(CERNLIB_DOUBLE)
-      DOUBLE PRECISION DA,DB,DC,DVAL
-#endif
-C          Convert to double precision
-      DA=A
-      DB=B
-      DC=C
-      DVAL=(DA-(DB+DC))*(DA+(DB+DC))*(DA-(DB-DC))*(DA+(DB-DC))
-C          Convert back to single precision
-      VAL=DVAL
-      DBLPCM=SQRT(ABS(VAL))/(2.*A)
-      RETURN
-      END
diff --git a/ISAJET/code/dblvec.F b/ISAJET/code/dblvec.F
deleted file mode 100644 (file)
index a2b657b..0000000
+++ /dev/null
@@ -1,28 +0,0 @@
-#include "isajet/pilot.h"
-      SUBROUTINE DBLVEC(P,DP)
-C
-C          Calculate double precision vector DP for 5-vector P.
-C          Exact components are 1,2,5 and larger of +,-
-C          Ver 6.44: Always use this, even if IF=SINGLE.
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      REAL P(5)
-      DOUBLE PRECISION DP(5),DPPL,DPMN
-      INTEGER K
-C
-      DO 100 K=1,5
-100   DP(K)=P(K)
-      IF(DP(4)+ABS(DP(3)).EQ.0.) RETURN
-      IF(DP(3).GT.0.) THEN
-        DPPL=DP(4)+DP(3)
-        DPMN=(DP(1)**2+DP(2)**2+DP(5)**2)/DPPL
-      ELSE
-        DPMN=DP(4)-DP(3)
-        DPPL=(DP(1)**2+DP(2)**2+DP(5)**2)/DPMN
-      ENDIF
-      DP(3)=0.5D0*(DPPL-DPMN)
-      DP(4)=0.5D0*(DPPL+DPMN)
-      RETURN
-      END
diff --git a/ISAJET/code/dboost.F b/ISAJET/code/dboost.F
deleted file mode 100644 (file)
index 373b4ab..0000000
+++ /dev/null
@@ -1,32 +0,0 @@
-#include "isajet/pilot.h"
-      SUBROUTINE DBOOST(ISIGN,F,P)
-C
-C          DOUBLE PRECISION BOOST OF 5-VECTOR P BY 5-VECTOR F WITH SIGN
-C          OF ISIGN. EXACT COMPONENTS ARE 1,2,5 AND LARGER OF +,-
-C
-      DIMENSION F(5),P(5)
-      DOUBLE PRECISION DF(5),DFPL,DFMN,DP(5),DPPL,DPMN,DBP,DSIGN
-C          COPY TO DOUBLE PRECISION
-      DO 100 K=1,5
-      DF(K)=F(K)
-100   DP(K)=P(K)
-      IF(ISIGN.GT.0) THEN
-        DSIGN=1.D0
-      ELSE
-        DSIGN=-1.D0
-      ENDIF
-C          PUT ON DOUBLE PRECISION SHELL
-      CALL DBLVEC(P,DP)
-C          BOOST
-      DBP=0.D0
-      DO 110 K=1,3
-110   DBP=DBP+DF(K)*DP(K)
-      DBP=DBP/DF(5)
-      DO 120 K=1,3
-120   DP(K)=DP(K)+DSIGN*DF(K)*DP(4)/DF(5)+DF(K)*DBP/(DF(4)+DF(5))
-      DP(4)=DF(4)*DP(4)/DF(5)+DSIGN*DBP
-C          COPY BACK
-      DO 130 K=1,4
-130   P(K)=DP(K)
-      RETURN
-      END
diff --git a/ISAJET/code/decay.F b/ISAJET/code/decay.F
deleted file mode 100644 (file)
index f03d2dd..0000000
+++ /dev/null
@@ -1,293 +0,0 @@
-#include "isajet/pilot.h"
-      SUBROUTINE DECAY(IP)
-C
-C          Decay particle IP from /PARTCL/ using /DKYTAB/ branching
-C          ratios and add decay products to /PARTCL/ with IORIG=IP.
-C          Forced decay modes are flagged by LOOK<0.
-C
-C          Auxiliary routines:
-C          DECPS1: generate masses for phase space
-C          DECPS2: generate 2-body decays and boosts for phase space
-C          DECVA:  V-A matrix elements
-C          DECTAU: tau decay matrix elements with polarization
-C          DECSS3: 3-body SUSY matrix element using /DKYSS3/
-C          DECJET: Hadronize partons from decay.
-C
-C          Matrix element for Dalitz decays and W mass for TP -> W BT
-C          are generated explicitly. W width is included.
-C
-C          Requirements for decay modes:
-C          (1) For Dalitz decays, particle 1 must be GM.
-C          (2) For V-A quark or lepton decays, particles 1 and 2 must
-C              be from (virtual) W.
-C          (3) For any decay into quarks, they must appear last.
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-#include "isajet/itapes.inc"
-#include "isajet/wcon.inc"
-#include "isajet/partcl.inc"
-#include "isajet/dkytab.inc"
-#include "isajet/jetset.inc"
-#include "isajet/jwork.inc"
-#include "isajet/const.inc"
-#include "isajet/primar.inc"
-#include "isajet/idrun.inc"
-#include "isajet/force.inc"
-#include "isajet/sstype.inc"
-#include "isajet/dkyss3.inc"
-C
-      REAL PGEN(5,5),BETA(3),REDUCE(5),WPROP,Z,TRY,RANF,AMASS,TWOME
-      REAL PSUM(5),SUM,PREST(4,6),DOT,PCM
-      REAL AMEE,REE,WTEE,SWAP,WT,A,B,C,GAMMA
-      REAL SMAX,SMIN,SVAL,TANMAX,TANMIN,TANVAL
-      LOGICAL WDECAY,DECVA,DECTAU,DECJET
-      INTEGER IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX,IPOINT,ID1,I1,I2
-      INTEGER NADD,NSTART,NEW,NADD1,J,IP,I,IDABS(5)
-      INTEGER K,JETIP,IDANTI,NPASS,MEIP,MEA
-      REAL DBLPCM,DECSS3,VAL
-C
-      DATA REDUCE/1.,1.,2.,5.,15./
-      DATA PSUM/5*0./
-      DATA TWOME/1.022006E-3/
-      DATA PREST/24*0./
-C
-C          Function definitions.
-C          Use double precision for PCM on 32-bit machines
-C
-#if defined(CERNLIB_SINGLE)
-      PCM(A,B,C)=SQRT((A**2-B**2-C**2)**2-(2.*B*C)**2)/(2.*A)
-#endif
-#if defined(CERNLIB_DOUBLE)
-      PCM(A,B,C)=DBLPCM(A,B,C)
-#endif
-      DOT(I1,I2)=PREST(4,I1)*PREST(4,I2)-PREST(1,I1)*PREST(1,I2)
-     $-PREST(2,I1)*PREST(2,I2)-PREST(3,I1)*PREST(3,I2)
-C          Charged W propagator.
-      WPROP(Z)=(Z-WMASS(2)**2)**2+(WMASS(2)*WGAM(2))**2
-C----------------------------------------------------------------------
-C          Select decay mode. Note IDENT(NPTCL+1)...IDENT(NPTCL+5)
-C          are always defined even if zero.
-C----------------------------------------------------------------------
-      IF(IDCAY(IP).NE.0) RETURN
-      IDLV1=IDENT(IP)
-      CALL FLAVOR(IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX)
-C          FLAVOR returns 0 for quark, but want IFL3=6 for top
-      IF(IABS(IDLV1).LT.10) IFL3=IDLV1
-      NPASS=0
-1     CONTINUE
-      NPASS=NPASS+1
-      WDECAY=.FALSE.
-      IF(NPASS.GT.NTRIES) GO TO 9998
-      IPOINT=LOOK(INDEX)
-      IF(IPOINT.EQ.0) RETURN
-C          IPOINT<0 flags a forced decay.
-      IF(IPOINT.LT.0) THEN
-        I=1
-        IF(IDENT(IP).LT.0) I=2
-        IPOINT=LOOK2(I,IABS(IPOINT))
-      ENDIF
-C
-C          Select mode.
-C
-      TRY=RANF()
-      IPOINT=IPOINT-1
-100   IPOINT=IPOINT+1
-      IF(TRY.GT.CBR(IPOINT)) GO TO 100
-      NADD=0
-      SUM=0.
-      NSTART=NPTCL+1
-      IF(NPTCL+5.GT.MXPTCL) GO TO 9999
-C
-C          Set up masses and IDENT codes.
-C
-      MEIP=MELEM(IPOINT)
-      DO 110 I=1,5
-        NEW=NPTCL+I
-        IDENT(NEW)=MODE(I,IPOINT)
-        IDABS(I)=IABS(IDENT(NEW))
-        IF(MODE(I,IPOINT).EQ.0) GO TO 110
-        NADD=NADD+1
-        IDLV1=IDENT(NEW)
-        PPTCL(5,NEW)=AMASS(IDLV1)
-        SUM=SUM+PPTCL(5,NEW)
-110   CONTINUE
-      NADD1=NADD-1
-      DO 120 J=1,5
-        PGEN(J,1)=PPTCL(J,IP)
-120   CONTINUE
-      PGEN(5,NADD)=PPTCL(5,NPTCL+NADD)
-C----------------------------------------------------------------------
-C          Carry out appropriate decay
-C----------------------------------------------------------------------
-C
-C          1-body decays.
-C
-      IF(NADD.EQ.1) THEN
-        DO 200 J=1,5
-          PPTCL(J,NPTCL+1)=PPTCL(J,IP)
-200     CONTINUE
-        GO TO 300
-      ENDIF
-C
-C          2-body phase space decays
-C
-      IF(NADD.EQ.2.AND.MEIP.EQ.0) THEN
-        CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
-        GO TO 300
-      ENDIF
-C
-C          N-body phase space decays
-C
-      IF(NADD.GT.2.AND.MEIP.EQ.0) THEN
-        CALL DECPS1(IP,NADD,PGEN)
-        CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
-        GO TO 300
-      ENDIF
-C
-C          Dalitz decays
-C
-      IF(NADD.EQ.3.AND.MEIP.EQ.1) THEN
-210     AMEE=TWOME*(PPTCL(5,IP)/TWOME)**RANF()
-        REE=(TWOME/AMEE)**2
-        WTEE=(1.-(AMEE/PPTCL(5,IP))**2)**3*SQRT(1.-REE)*(1.+.5*REE)
-        IF(WTEE.LT.RANF()) GO TO 210
-        PGEN(5,2)=AMEE
-        CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
-        GO TO 300
-      ENDIF
-C
-C          omega/phi decays (for reasons lost in history...)
-C
-      IF(NADD.EQ.3.AND.MEIP.EQ.2) THEN
-220     CALL DECPS1(IP,NADD,PGEN)
-        CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
-        WT=(PPTCL(5,NPTCL+1)*PPTCL(5,NPTCL+2)*PPTCL(5,NPTCL+3))**2
-     $  -(PPTCL(5,NPTCL+1)*DOT(2,3))**2
-     $  -(PPTCL(5,NPTCL+2)*DOT(1,3))**2
-     $  -(PPTCL(5,NPTCL+3)*DOT(1,2))**2
-     $  +2.*DOT(1,2)*DOT(2,3)*DOT(1,3)
-        IF(WT.LT.RANF()*PPTCL(5,IP)**6/108.) GO TO 220
-        GO TO 300
-      ENDIF
-C
-C          V-A decays
-C
-      IF(NADD.EQ.3.AND.MEIP.EQ.3) THEN
-230     CALL DECPS1(IP,NADD,PGEN)
-        CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
-        IF(.NOT.DECVA(IP,NADD,IDABS,PREST)) GO TO 230
-        GO TO 300
-      ENDIF
-C
-C          Top decays
-C          Generate mass for TP -> W BT with Breit-Wigner. 
-C          W couples to 1+2 so swap 1<->3. Then m2+m3 < m < m0-m1.
-C
-      IF(NADD.EQ.3.AND.MEIP.EQ.4) THEN
-        WDECAY=.TRUE.
-        SWAP=PPTCL(5,NPTCL+1)
-        PPTCL(5,NPTCL+1)=PPTCL(5,NPTCL+3)
-        PPTCL(5,NPTCL+3)=SWAP
-        SMAX=(PPTCL(5,IP)-PPTCL(5,NPTCL+1))**2
-        SMIN=(PPTCL(5,NPTCL+2)+PPTCL(5,NPTCL+3))**2
-        TANMAX=ATAN((SMAX-WMASS(2)**2)/(WMASS(2)*WGAM(2)))
-        TANMIN=ATAN((SMIN-WMASS(2)**2)/(WMASS(2)*WGAM(2)))
-240     TANVAL=RANF()*(TANMAX-TANMIN)+TANMIN
-        SVAL=WMASS(2)**2+WMASS(2)*WGAM(2)*TAN(TANVAL)
-        IF(SVAL.LT.SMIN.OR.SVAL.GT.SMAX) GO TO 240
-        PGEN(5,2)=SQRT(SVAL)
-        PGEN(5,3)=PPTCL(5,NPTCL+3)
-        CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
-        IF(.NOT.DECVA(IP,NADD,IDABS,PREST)) GO TO 240
-        DO 241 K=1,5
-          SWAP=PPTCL(K,NPTCL+1)
-          PPTCL(K,NPTCL+1)=PPTCL(K,NPTCL+3)
-          PPTCL(K,NPTCL+3)=SWAP
-241     CONTINUE
-        PGEN(5,3)=PPTCL(5,NPTCL+3)
-        DO 242 K=1,4
-          SWAP=PREST(K,1)
-          PREST(K,1)=PREST(K,3)
-          PREST(K,3)=SWAP
-242     CONTINUE
-        GO TO 300
-      ENDIF
-C
-C          TAU decays. These are special because they take polarization
-C          into account.
-C
-      IF(MEIP.EQ.5.OR.MEIP.EQ.6.OR.MEIP.EQ.7) THEN
-250     CALL DECPS1(IP,NADD,PGEN)
-        CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
-        IF(.NOT.DECTAU(IP,NADD,MEIP,IDABS,PREST)) GO TO 250
-        GO TO 300
-      ENDIF
-C
-C          3-body SUSY decays
-C
-      IF(MEIP.LT.0.AND.NADD.EQ.3) THEN
-        MEA=IABS(MEIP)
-        IF(WTSS3(MEA).LE.0) THEN
-          DO 260 I=1,1000
-            CALL DECPS1(IP,NADD,PGEN)
-            CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
-            VAL=DECSS3(IP,MEA)
-            WTSS3(MEA)=MAX(WTSS3(MEA),VAL)
-260       CONTINUE
-          IF(WTSS3(MEA).LE.0) GO TO 9998
-        ENDIF
-261     CALL DECPS1(IP,NADD,PGEN)
-        CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
-        VAL=DECSS3(IP,MEA)
-        WTSS3(MEA)=MAX(WTSS3(MEA),VAL)
-        IF(VAL.LT.WTSS3(MEA)*RANF()) GO TO 261
-        GO TO 300
-      ENDIF
-C
-C          Should never fall through
-C
-      GO TO 9998
-C----------------------------------------------------------------------
-C          Swap particles and antiparticles if IDENT(IP)<0
-C          Note forced modes for antiparticles are conjugated in table.
-C----------------------------------------------------------------------
-300   CONTINUE
-      IF(IDENT(IP).LT.0.AND.IDENT(IP).NE.-20) THEN
-        DO 310 I=1,NADD
-          ID1=IDENT(NPTCL+I)
-          IDENT(NPTCL+I)=IDANTI(ID1)
-310     CONTINUE
-      ENDIF
-C
-C          Set IORIG and IDCAY.
-C
-      NPTCL=NPTCL+NADD
-      IDCAY(IP)=IPACK*NSTART+NPTCL
-      JETIP=IABS(IORIG(IP))/IPACK
-      DO 320 I=NSTART,NPTCL
-        IORIG(I)=IP
-        IDCAY(I)=0
-320   CONTINUE
-C
-C          Evolve and hadronize partons. If it fails, start over.
-C
-      IF(IDABS(NADD).LT.10.OR.MOD(IDENT(NPTCL),100).EQ.0) THEN
-        IF(.NOT.DECJET(IP,NADD,IDABS,PREST,WDECAY,BETA,GAMMA))
-     $  GO TO 1
-      ENDIF
-C
-      RETURN
-C----------------------------------------------------------------------
-C          Error messages.
-C----------------------------------------------------------------------
-9999  CALL PRTEVT(0)
-      WRITE(ITLIS,99990) NPTCL
-99990 FORMAT(//5X,'ERROR IN DECAY...NPTCL > ',I6)
-      RETURN
-9998  CALL PRTEVT(0)
-      WRITE(ITLIS,99980) IP
-99980 FORMAT(//5X,'ERROR IN DECAY...NO DECAY FOUND FOR PARTICLE',I6)
-      RETURN
-      END
diff --git a/ISAJET/code/decjet.F b/ISAJET/code/decjet.F
deleted file mode 100644 (file)
index 03454b1..0000000
+++ /dev/null
@@ -1,380 +0,0 @@
-#include "isajet/pilot.h"
-      LOGICAL FUNCTION DECJET(IP,NADD,IDABS,PREST,WDECAY,BETA,GAMMA)
-C
-C          Auxiliary routine for DECAY. Evolve and hadronize partons.
-C          Check conservation laws. Return TRUE if OK, FALSE otherwise.
-C
-C          IP = particle to be decayed.
-C          NADD = number of products (NPTCL+1, ..., NPTCL+NADD).
-C          IDABS = absolute values of decay IDENT's.
-C          PREST = 4-momenta in rest frame.
-C          WDECAY = logical flag for real W decay.
-C          BETA,GAMMA = boost parameters.
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-#include "isajet/itapes.inc"
-#include "isajet/wcon.inc"
-#include "isajet/partcl.inc"
-#include "isajet/dkytab.inc"
-#include "isajet/jetset.inc"
-#include "isajet/jwork.inc"
-#include "isajet/const.inc"
-C   
-      REAL PGEN(5,5),RND(5),U(3),BETA(3),IDQK(3),ROT(3,3),PSAVE(3)  
-     1,REDUCE(5),WPROP,Z,TRY,RANF,AMASS,TWOME,CHARGE    
-      REAL PSUM(5),POLD(4),PNEW(4),SUM,WTMAX,SUM1,SUM2  
-      REAL PREST(4,6),PWREST(5),BETAW(3),DOT,PCM    
-      REAL AMEE,REE,WTEE,SWAP,RNEW,WT,QCM,PHI,S12,S12MAX,GAMMAW,BP  
-      REAL PJET,CTHQK,STHQK,CPHIQK,SPHIQK,SUMQ,A,B,C,GAMMA  
-      REAL CHARGW   
-      LOGICAL WDECAY    
-      INTEGER IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX,IPOINT,ID1,I1,I2,I3  
-      INTEGER NADD,NSTART,NEW,NADD1,J,IP,I,IDABS(5),NEXT    
-      INTEGER JJ1,II,K1,K,NJSAVE,NJSAV1,NJSAV2,NJ1,NPRTN,KK,NHDRN1  
-      INTEGER IFAIL,JSAVE,JETIP,JET,NJADD,NPTLV1,IDANTI,NPJET(5)    
-      INTEGER NHDRN,NPJET3,NPTCLW,NPBEG(5)  
-C   
-C          Copy decay products into /JETSET/ and do QCD evolution.  
-C   
-      IF(NJSET+NADD.GT.MXJSET) GO TO 9998   
-      NJSAVE=NJSET  
-      NSTART=NPTCL-NADD+1   
-      NPTCL=NSTART-1    
-      DO 100 I=1,NADD   
-        NJSET=NJSET+1   
-        DO 110 K=1,4    
-110     PJSET(K,NJSET)=PREST(K,I)   
-        PJSET(5,NJSET)=PPTCL(5,NPTCL+I) 
-        JORIG(NJSET)=JPACK*I    
-        JTYPE(NJSET)=IDENT(NPTCL+I) 
-        JDCAY(NJSET)=0  
-        JMATCH(NJSET)=JPACK*(NJSAVE+1)+NJSAVE+NADD  
-100   CONTINUE  
-C   
-C          For heavy quarks match 1+2 and 3+(1+2). Boost 1+2 to rest.   
-C   
-      IF(WDECAY) THEN   
-        JMATCH(NJSAVE+1)=NJSAVE+2   
-        JMATCH(NJSAVE+2)=NJSAVE+1   
-        NJSET=NJSET+1   
-        DO 120 K=1,4    
-          PWREST(K)=PJSET(K,NJSAVE+1)+PJSET(K,NJSAVE+2) 
-          PJSET(K,NJSET)=PWREST(K)  
-120     CONTINUE    
-        PWREST(5)=SQRT(PWREST(4)**2-PWREST(1)**2-PWREST(2)**2   
-     $  -PWREST(3)**2)  
-        PJSET(5,NJSET)=PWREST(5)    
-        JMATCH(NJSAVE+3)=NJSAVE+4   
-        JMATCH(NJSAVE+4)=NJSAVE+3   
-        JORIG(NJSAVE+4)=-1  
-        IDLV1=JTYPE(NJSAVE+1)   
-        CHARGW=CHARGE(IDLV1)    
-        IDLV1=JTYPE(NJSAVE+2)   
-        CHARGW=CHARGW+CHARGE(IDLV1) 
-        JTYPE(NJSAVE+4)=80*SIGN(1.,CHARGW)  
-        JDCAY(NJSAVE+4)=0   
-C          Boost W vectors to rest. 
-        DO 130 K=1,3    
-130     BETAW(K)=PWREST(K)/PWREST(4)    
-        GAMMAW=PWREST(4)/PWREST(5)  
-        NJSAV1=NJSAVE+1 
-        NJSAV2=NJSAVE+2 
-        DO 140 J=NJSAV1,NJSAV2  
-          BP=BETAW(1)*PJSET(1,J)+BETAW(2)*PJSET(2,J)+BETAW(3)*PJSET(3,J)    
-          DO 141 K=1,3  
-141       PJSET(K,J)=PJSET(K,J)-GAMMAW*BETAW(K)*(PJSET(4,J) 
-     $    -BP*GAMMAW/(GAMMAW+1.))   
-          PJSET(4,J)=GAMMAW*(PJSET(4,J)-BP) 
-140     CONTINUE    
-      ENDIF 
-C   
-C          Do evolution and save new W momentum. Start from parent  
-C          mass or NADD*energy. 
-      NJSAV1=NJSAVE+1   
-      DO 150 J=NJSAV1,NJSET 
-        IF(IABS(JTYPE(J)).LT.10.OR.MOD(JTYPE(J),100).EQ.0) THEN 
-          JDCAY(J)=-1   
-          PJSET(5,J)=AMIN1(PPTCL(5,IP),NADD*PJSET(4,J)) 
-        ENDIF   
-150   CONTINUE  
-C   
-      CALL QCDJET(NJSAVE+1) 
-C   
-      IF(WDECAY) THEN   
-        PWREST(4)=PJSET(4,NJSAVE+4) 
-        GAMMAW=PWREST(4)/PWREST(5)  
-        DO 200 K=1,3    
-          PWREST(K)=PJSET(K,NJSAVE+4)   
-          BETAW(K)=PWREST(K)/PWREST(4)  
-200     CONTINUE    
-      ENDIF 
-C   
-C          Put final partons in particle table - temporary IORIG.   
-C          Also include virtual or real W momentum for quark decays.    
-C   
-      NJ1=NJSAVE+1  
-      IF(WDECAY) THEN   
-C          Real or virtual W.   
-        NPTCL=NPTCL+1   
-        NPTCLW=NPTCL    
-        DO 210 K=1,5    
-210     PPTCL(K,NPTCL)=PJSET(K,NJSAVE+4)    
-        IORIG(NPTCL)=IP 
-        IDENT(NPTCL)=JTYPE(NJSAVE+4)    
-        IDCAY(NPTCL)=0  
-C          Jet 3    
-        NPBEG(3)=NPTCL+1    
-        DO 220 J=NJ1,NJSET  
-          IF(JDCAY(J).NE.0) GO TO 220   
-          IF(JORIG(J)/JPACK.NE.3) GO TO 220 
-          NPTCL=NPTCL+1 
-          DO 221 K=1,5  
-221       PPTCL(K,NPTCL)=PJSET(K,J) 
-          IORIG(NPTCL)=3*IPACK+IP   
-          IDENT(NPTCL)=JTYPE(J) 
-          IDCAY(NPTCL)=0    
-220     CONTINUE    
-C          Jets 1 and 2 
-        NPJET3=NPTCL    
-        DO 230 JET=1,2  
-          NPBEG(JET)=NPTCL+1    
-          DO 240 J=NJ1,NJSET    
-            IF(JDCAY(J).NE.0) GO TO 240 
-            IF(JORIG(J)/JPACK.NE.JET) GO TO 240 
-            NPTCL=NPTCL+1   
-            BP=BETAW(1)*PJSET(1,J)+BETAW(2)*PJSET(2,J)  
-     $      +BETAW(3)*PJSET(3,J)    
-            DO 241 K=1,3    
-241         PPTCL(K,NPTCL)=PJSET(K,J)+GAMMAW*BETAW(K)*(PJSET(4,J)   
-     $      +BP*GAMMAW/(GAMMAW+1.)) 
-            PPTCL(4,NPTCL)=GAMMAW*(PJSET(4,J)+BP)   
-            PPTCL(5,NPTCL)=PJSET(5,J)   
-            IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+NPTCLW  
-            IDENT(NPTCL)=JTYPE(J)   
-            IDCAY(NPTCL)=0  
-240       CONTINUE  
-230     CONTINUE    
-C          Quark decays to W plus jet 3; then W decays. 
-        IDCAY(IP)=IPACK*NPTCLW+NPJET3   
-        IDCAY(NPTCLW)=IPACK*(NPJET3+1)+NPTCL    
-      ELSE  
-C          Not quark decay, so just copy partons.   
-        DO 250 JET=1,NADD   
-          NPBEG(JET)=NPTCL+1    
-          DO 260 J=NJ1,NJSET    
-            IF(JDCAY(J).NE.0) GO TO 260 
-            IF(JORIG(J)/JPACK.NE.JET) GO TO 260 
-            NPTCL=NPTCL+1   
-            DO 261 K=1,5    
-261         PPTCL(K,NPTCL)=PJSET(K,J)   
-            IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+IP  
-            IDENT(NPTCL)=JTYPE(J)   
-            IDCAY(NPTCL)=0  
-260       CONTINUE  
-250     CONTINUE    
-        IDCAY(IP)=NSTART*IPACK+NPTCL    
-      ENDIF 
-      NHDRN=NPTCL   
-C   
-C          Hadronize quarks and rotate to proper angles.    
-C   
-      DO 300 JET=1,NADD 
-        NPRTN=NPBEG(JET)-1  
-        DO 310 I=NJ1,NJSET  
-          IF(JDCAY(I).NE.0) GO TO 310   
-          IF(JORIG(I)/JPACK.NE.JET) GO TO 310   
-          NPRTN=NPRTN+1 
-          IF(IABS(JTYPE(I)).GE.10.AND.MOD(JTYPE(I),100).NE.0)   
-     $    GO TO 330 
-C   
-C          Fragment parton: 
-          NEXT=NPTCL+1  
-          PJET=SQRT(PJSET(1,I)**2+PJSET(2,I)**2+PJSET(3,I)**2)  
-          CTHQK=PJSET(3,I)/PJET 
-          STHQK=1.-CTHQK**2
-          IF(STHQK.LT.1) THEN
-            STHQK=SQRT(STHQK)
-            CPHIQK=PJSET(1,I)/(PJET*STHQK)
-            SPHIQK=PJSET(2,I)/(PJET*STHQK)
-          ELSE
-            STHQK=0
-            CPHIQK=1
-            SPHIQK=0
-          ENDIF
-          CALL JETGEN(I)    
-          IF(NEXT.GT.NPTCL) GO TO 310   
-          ROT(1,1)=CPHIQK*CTHQK 
-          ROT(2,1)=SPHIQK*CTHQK 
-          ROT(3,1)=-STHQK   
-          ROT(1,2)=-SPHIQK  
-          ROT(2,2)=CPHIQK   
-          ROT(3,2)=0.   
-          ROT(1,3)=CPHIQK*STHQK 
-          ROT(2,3)=SPHIQK*STHQK 
-          ROT(3,3)=CTHQK    
-C   
-          DO 320 II=NEXT,NPTCL  
-            DO 321 K=1,3    
-              PSAVE(K)=PPTCL(K,II)  
-              PPTCL(K,II)=0.    
-321         CONTINUE    
-            DO 322 K=1,3    
-            DO 322 KK=1,3   
-322         PPTCL(K,II)=PPTCL(K,II)+ROT(K,KK)*PSAVE(KK) 
-            IORIG(II)=IPACK*JET+NPRTN   
-            IDCAY(II)=0 
-320       CONTINUE  
-          IDCAY(NPRTN)=NEXT*IPACK+NPTCL 
-          GO TO 310 
-C   
-C          or add lepton:   
-330       NPTCL=NPTCL+1 
-          DO 331 K=1,5  
-331       PPTCL(K,NPTCL)=PJSET(K,I) 
-          IORIG(NPTCL)=IPACK*JET+NPRTN  
-          IDENT(NPTCL)=JTYPE(I) 
-          IDCAY(NPTCL)=0    
-          IDCAY(NPRTN)=NPTCL*IPACK+NPTCL    
-310     CONTINUE    
-        NPJET(JET)=NPTCL    
-300   CONTINUE  
-C   
-C          Reset NJSET so decay jets do not appear in /JETSET/  
-      NJADD=NJSET   
-      NJSET=NJSAVE  
-C   
-C          Check for at least two particles 
-      IF(NPTCL.LT.NHDRN+2) THEN 
-        NPTCL=NSTART-1  
-        DECJET=.FALSE.  
-        RETURN  
-      ENDIF 
-C   
-C          Conserve charge  
-C   
-      SUMQ=0.   
-      NHDRN1=NHDRN+1    
-      DO 400 I=NHDRN1,NPTCL 
-        IDLV1=IDENT(I)  
-        SUMQ=SUMQ+CHARGE(IDLV1) 
-400   CONTINUE  
-      IDLV1=IDENT(IP)   
-      SUMQ=SUMQ-CHARGE(IDLV1)   
-C   
-      IF(SUMQ.EQ.0.) GO TO 500  
-C   
-C          Charge wrong--fix it by swapping UP and DN quarks.   
-      DO 410 I=NHDRN1,NPTCL 
-        ID1=IDENT(I)    
-        IF(IABS(ID1).GT.1000) GO TO 410 
-        I1=MOD(IABS(ID1)/100,10)    
-        I2=MOD(IABS(ID1)/10,10) 
-        I3=MOD(IABS(ID1),10)    
-        IF(I1.EQ.1.AND.I2.GT.2.AND.SUMQ*ID1.GT.0.) THEN 
-          IDENT(I)=ISIGN(200+10*I2+I3,ID1)  
-        ELSEIF(I1.EQ.2.AND.I2.GT.2.AND.SUMQ*ID1.LT.0.) THEN 
-          IDENT(I)=ISIGN(100+10*I2+I3,ID1)  
-        ELSEIF(I1.EQ.1.AND.I2.EQ.2.AND.SUMQ*ID1.GT.0.) THEN 
-          IDENT(I)=110+I3   
-        ELSEIF(I1.EQ.1.AND.I2.EQ.1) THEN    
-          IDENT(I)=(120+I3)*(-SIGN(1.,SUMQ))    
-        ELSE    
-          GO TO 410 
-        ENDIF   
-        SUMQ=SIGN(ABS(SUMQ)-1.,SUMQ)    
-        IDLV1=IDENT(I)  
-        PPTCL(5,I)=AMASS(IDLV1) 
-        PPTCL(4,I)=SQRT(PPTCL(1,I)**2+PPTCL(2,I)**2+PPTCL(3,I)**2   
-     $  +PPTCL(5,I)**2) 
-C          Sum cannot vanish for fractionally charged initial particle. 
-        IF(ABS(SUMQ).LT.1.) GO TO 500   
-410   CONTINUE  
-C          Failed to conserve charge.   
-      NPTCL=NSTART-1    
-      DECJET=.FALSE.    
-      RETURN    
-C   
-C          Rescale momenta for correct mass 
-C   
-500   CONTINUE  
-      IF(WDECAY) THEN   
-C          Quark decay. First rescale jet3 + W  
-        DO 510 K=1,5    
-510     PPTCL(K,NPTCL+1)=PPTCL(K,NPTCLW)    
-        NPTLV1=NPTCL+1  
-        DO 520 K=1,3    
-520     PSUM(K)=0.  
-        PSUM(4)=PPTCL(5,IP) 
-        PSUM(5)=PSUM(4) 
-        CALL RESCAL(NPJET(2)+1,NPTLV1,PSUM,IFAIL) 
-        IF(IFAIL.NE.0) THEN 
-          NPTCL=NSTART-1    
-          DECJET=.FALSE.    
-          RETURN    
-        ENDIF   
-        DO 530 K=1,3    
-530     BETAW(K)=PPTCL(K,NPTCL+1)/PPTCL(4,NPTCL+1)  
-        GAMMAW=PPTCL(4,NPTCL+1)/PPTCL(5,NPTCL+1)    
-C          Then rescale W   
-        PSUM(4)=PPTCL(5,NPTCLW) 
-        PSUM(5)=PSUM(4) 
-        CALL RESCAL(NHDRN1,NPJET(2),PSUM,IFAIL) 
-        IF(IFAIL.NE.0) THEN 
-          NPTCL=NSTART-1    
-          DECJET=.FALSE.    
-          RETURN    
-        ENDIF   
-      ELSE  
-C          General decay with no W. 
-        DO 550 K=1,3    
-550     PSUM(K)=0.  
-        PSUM(4)=PPTCL(5,IP) 
-        PSUM(5)=PSUM(4) 
-        NPTLV1=NPTCL    
-        CALL RESCAL(NHDRN1,NPTLV1,PSUM,IFAIL)   
-        IF(IFAIL.NE.0) THEN 
-          NPTCL=NSTART-1    
-          DECJET=.FALSE.    
-          RETURN    
-        ENDIF   
-      ENDIF 
-C   
-C          Boost back to lab frame. Reset IORIG.    
-C   
-      IF(WDECAY) THEN   
-        DO 600 I=NHDRN1,NPTCL  
-          JET=IORIG(I)/IPACK    
-          IF(JET.NE.1.AND.JET.NE.2) GO TO 600   
-          BP=BETAW(1)*PPTCL(1,I)+BETAW(2)*PPTCL(2,I)+BETAW(3)*PPTCL(3,I)    
-          DO 610 J=1,3  
-610       PPTCL(J,I)=PPTCL(J,I)+GAMMAW*BETAW(J)*(PPTCL(4,I) 
-     $    +BP*GAMMAW/(GAMMAW+1.))   
-          PPTCL(4,I)=GAMMAW*(PPTCL(4,I)+BP) 
-600     CONTINUE    
-      ENDIF 
-C   
-      DO 620 I=NSTART,NPTCL 
-        IORIG(I)=MOD(IORIG(I),IPACK)    
-        BP=BETA(1)*PPTCL(1,I)+BETA(2)*PPTCL(2,I)+BETA(3)*PPTCL(3,I) 
-        DO 621 J=1,3    
-          PPTCL(J,I)=PPTCL(J,I)+GAMMA*BETA(J)*(PPTCL(4,I)   
-     $    +BP*GAMMA/(GAMMA+1.)) 
-621     CONTINUE    
-        PPTCL(4,I)=GAMMA*(PPTCL(4,I)+BP)    
-620   CONTINUE  
-C   
-C          Normal exit  
-C   
-      DECJET=.TRUE. 
-      RETURN    
-C   
-C          Error messages.  
-C   
-9998  DECJET=.FALSE.
-      CALL PRTEVT(0)    
-      WRITE(ITLIS,99980) NJSET  
-99980 FORMAT(//5X,'ERROR IN DECJET...NJSET > ',I5)  
-      RETURN    
-      END   
diff --git a/ISAJET/code/decps1.F b/ISAJET/code/decps1.F
deleted file mode 100644 (file)
index 1b5131e..0000000
+++ /dev/null
@@ -1,75 +0,0 @@
-#include "isajet/pilot.h"
-      SUBROUTINE DECPS1(IP,NADD,PGEN)
-C
-C          Generate masses for uniform NADD-body phase space in DECPS2.
-C          Auxiliary routine for DECAY.
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-C
-#include "isajet/itapes.inc"
-#include "isajet/partcl.inc"
-C
-      INTEGER IP,NADD
-      REAL PGEN(5,5)
-      REAL REDUCE(5),RND(5)
-      REAL RANF,PCM,DBLPCM
-      REAL WTMAX,SUM1,SUM2,SUM,RNEW,WT,A,B,C
-      INTEGER I,NADD1,J,I1,JJ1,JSAVE
-C
-C          Function definitions.
-C
-#if defined(CERNLIB_SINGLE)
-      PCM(A,B,C)=SQRT((A-B-C)*(A+B+C)*(A-B+C)*(A+B-C))/(2.*A)
-#endif
-#if defined(CERNLIB_DOUBLE)
-      PCM(A,B,C)=DBLPCM(A,B,C)
-#endif
-C
-      DATA REDUCE/1.,1.,2.,5.,15./
-C
-C          Calculate maximum phase-space weight.
-C
-      IF(NADD.LE.2) RETURN
-      NADD1=NADD-1
-      WTMAX=1./REDUCE(NADD)
-      SUM=0
-      DO 100 I=1,NADD
-        SUM=SUM+PPTCL(5,NPTCL+I)
-100   CONTINUE
-      SUM1=PGEN(5,1)
-      SUM2=SUM-PPTCL(5,NPTCL+1)
-      DO 110 I=1,NADD1
-        WTMAX=WTMAX*PCM(SUM1,SUM2,PPTCL(5,NPTCL+I))
-        SUM1=SUM1-PPTCL(5,NPTCL+I)
-        SUM2=SUM2-PPTCL(5,NPTCL+I+1)
-110   CONTINUE
-C
-C          Generate masses for uniform NADD-body phase space.
-C
-200   CONTINUE
-      RND(1)=1.
-      DO 210 I=2,NADD1
-        RNEW=RANF()
-        I1=I-1
-        DO 220 JJ1=1,I1
-          J=I-JJ1
-          JSAVE=J+1
-          IF(RNEW.LE.RND(J)) GO TO 210
-          RND(JSAVE)=RND(J)
-220     CONTINUE
-210   RND(JSAVE)=RNEW
-      RND(NADD)=0.
-      WT=1.
-      SUM1=SUM
-      DO 230 I=2,NADD
-        SUM1=SUM1-PPTCL(5,NPTCL+I-1)
-        PGEN(5,I)=SUM1+RND(I)*(PGEN(5,1)-SUM)
-        IF(PGEN(5,I-1).LE.PGEN(5,I)+PPTCL(5,NPTCL+I-1)) GO TO 200
-        WT=WT*PCM(PGEN(5,I-1),PGEN(5,I),PPTCL(5,NPTCL+I-1))
-230   CONTINUE
-      IF(WT.LT.RANF()*WTMAX) GO TO 200
-C
-      RETURN
-      END
diff --git a/ISAJET/code/decps2.F b/ISAJET/code/decps2.F
deleted file mode 100644 (file)
index 8c1a60d..0000000
+++ /dev/null
@@ -1,76 +0,0 @@
-#include "isajet/pilot.h"
-      SUBROUTINE DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
-C
-C          Carry out decays using masses from DECPS1 or special matrix
-C          elements.
-C          Auxiliary routine for DECAY.
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-C
-#include "isajet/itapes.inc"
-#include "isajet/partcl.inc"
-#include "isajet/const.inc"
-C
-      INTEGER IP,NADD
-      REAL PGEN(5,5),PREST(4,6)
-      REAL PCM,DBLPCM,RANF
-      REAL U(3),BETA(3)
-      REAL QCM,PHI,A,B,C,GAMMA,BP
-      INTEGER I,J,NADD1,II,K,K1
-C
-C          Function definitions.
-C
-#if defined(CERNLIB_SINGLE)
-      PCM(A,B,C)=SQRT((A-B-C)*(A+B+C)*(A-B+C)*(A+B-C))/(2.*A)
-#endif
-#if defined(CERNLIB_DOUBLE)
-      PCM(A,B,C)=DBLPCM(A,B,C)
-#endif
-C
-C          Carry out two-body decays in PGEN frames.
-C
-      NADD1=NADD-1
-100   CONTINUE
-      DO 110 I=1,NADD1
-        QCM=PCM(PGEN(5,I),PGEN(5,I+1),PPTCL(5,NPTCL+I))
-        U(3)=2.*RANF()-1.
-        PHI=2.*PI*RANF()
-        U(1)=SQRT(1.-U(3)**2)*COS(PHI)
-        U(2)=SQRT(1.-U(3)**2)*SIN(PHI)
-        DO 120 J=1,3
-          PPTCL(J,NPTCL+I)=QCM*U(J)
-          PGEN(J,I+1)=-PPTCL(J,NPTCL+I)
-120     CONTINUE
-        PPTCL(4,NPTCL+I)=SQRT(QCM**2+PPTCL(5,NPTCL+I)**2)
-        PGEN(4,I+1)=SQRT(QCM**2+PGEN(5,I+1)**2)
-110   CONTINUE
-C
-      DO 130 J=1,4
-        PPTCL(J,NPTCL+NADD)=PGEN(J,NADD)
-130   CONTINUE
-C
-C          Boost PGEN frames to lab frame, saving momenta in rest frame.
-C
-      DO 200 II=1,NADD1
-        I=NADD-II
-        DO 210 J=1,3
-          BETA(J)=PGEN(J,I)/PGEN(4,I)
-210     CONTINUE
-        GAMMA=PGEN(4,I)/PGEN(5,I)
-        DO 220 K=I,NADD
-          K1=NPTCL+K
-          BP=BETA(1)*PPTCL(1,K1)+BETA(2)*PPTCL(2,K1)+BETA(3)*PPTCL(3,K1)
-          DO 230 J=1,3
-            PREST(J,K)=PPTCL(J,K1)
-            PPTCL(J,K1)=PPTCL(J,K1)+GAMMA*BETA(J)*(PPTCL(4,K1)
-     $      +BP*GAMMA/(GAMMA+1.))
-230       CONTINUE
-          PREST(4,K)=PPTCL(4,K1)
-          PPTCL(4,K1)=GAMMA*(PPTCL(4,K1)+BP)
-220     CONTINUE
-200   CONTINUE
-C
-      RETURN
-      END
diff --git a/ISAJET/code/decss3.F b/ISAJET/code/decss3.F
deleted file mode 100644 (file)
index 5ab808f..0000000
+++ /dev/null
@@ -1,163 +0,0 @@
-#include "isajet/pilot.h"
-      FUNCTION DECSS3(IP,MEA)
-C
-C          Compute matrix element for mode MEA of particle IP using
-C          poles and couplings in /DKYSS3/.
-C          Auxiliary routine for DECAY.
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-C
-#include "isajet/itapes.inc"
-#include "isajet/partcl.inc"
-#include "isajet/const.inc"
-#include "isajet/dkyss3.inc"
-C
-      LOGICAL KIN(4),KINP(4)
-      INTEGER IP,MEA,I,J,JP,II,PTYPE1,PTYPE2
-      REAL DECSS3
-      REAL AM0SQ,AM1SQ,AM2SQ,AM3SQ,S12,S13,S23
-      REAL D12,D13,D23,D01,D02,D03,AS,BS,CS,DS,MSQ
-      REAL DOT4
-      COMPLEX A,B,C,D,AC,BC,CC,DC,AP,BP,CP,DP,APC,BPC,CPC,DPC,MMPD
-C
-      DOT4(I,J)=PPTCL(4,I)*PPTCL(4,J)-PPTCL(1,I)*PPTCL(1,J)-
-     $PPTCL(2,I)*PPTCL(2,J)-PPTCL(3,I)*PPTCL(3,J)
-C
-C          Kinematics
-C
-      AM0SQ=PPTCL(5,IP)**2
-      AM1SQ=PPTCL(5,NPTCL+1)**2
-      AM2SQ=PPTCL(5,NPTCL+2)**2
-      AM3SQ=PPTCL(5,NPTCL+3)**2
-      D12=DOT4(NPTCL+1,NPTCL+2)
-      D13=DOT4(NPTCL+1,NPTCL+3)
-      D23=DOT4(NPTCL+2,NPTCL+3)
-      D01=DOT4(IP,NPTCL+1)
-      D02=DOT4(IP,NPTCL+2)
-      D03=DOT4(IP,NPTCL+3)
-      S12=2*D12+AM1SQ+AM2SQ
-      S13=2*D13+AM1SQ+AM3SQ
-      S23=2*D23+AM2SQ+AM3SQ
-C
-C          Generic matrix element
-C
-C          Loop over diagrams
-      DECSS3=0.
-      DO J=J1SS3(MEA),J2SS3(MEA)
-       PTYPE1=KSS3(J)
-       A=ZISS3(1,J)
-       B=ZISS3(2,J)
-       C=ZFSS3(1,J)
-       D=ZFSS3(2,J)
-       AC=CONJG(A)
-       BC=CONJG(B)
-       CC=CONJG(C)
-       DC=CONJG(D)
-       AS=A*AC
-       BS=B*BC
-       CS=C*CC
-       DS=D*DC
-       DO JP=J,J2SS3(MEA)
-        MSQ=0.
-        DO II=1,4
-          KIN(II)=.FALSE.
-          KINP(II)=.FALSE.
-        END DO
-        IF ((PPTCL(5,IP)-PPTCL(5,NPTCL+1)).LT.AMSS3(J)) KIN(1)=.TRUE.
-        IF ((PPTCL(5,IP)-PPTCL(5,NPTCL+3)).LT.AMSS3(J)) KIN(2)=.TRUE.
-        IF ((PPTCL(5,IP)-PPTCL(5,NPTCL+2)).LT.AMSS3(J)) KIN(3)=.TRUE.
-        IF ((PPTCL(5,IP)-PPTCL(5,NPTCL+1)).LT.AMSS3(J)) KIN(4)=.TRUE.
-        IF ((PPTCL(5,IP)-PPTCL(5,NPTCL+1)).LT.AMSS3(JP)) KINP(1)=.TRUE.
-        IF ((PPTCL(5,IP)-PPTCL(5,NPTCL+3)).LT.AMSS3(JP)) KINP(2)=.TRUE.
-        IF ((PPTCL(5,IP)-PPTCL(5,NPTCL+2)).LT.AMSS3(JP)) KINP(3)=.TRUE.
-        IF ((PPTCL(5,IP)-PPTCL(5,NPTCL+1)).LT.AMSS3(JP)) KINP(4)=.TRUE.
-        IF (J.EQ.JP) THEN
-         IF (PTYPE1.EQ.1.AND.KIN(1)) THEN
-          MSQ=32*(((AS+BS)*(CS+DS)+4*REAL(A*BC*C*DC))*D03*D12+
-     $            ((AS+BS)*(CS+DS)-4*REAL(A*BC*C*DC))*D02*D13+
-     $             (BS-AS)*(CS+DS)*SQRT(AM0SQ*AM1SQ)*D23)/
-     $            (S23-AMSS3(J)**2)**2
-         ELSE IF (PTYPE1.EQ.2.AND.KIN(2)) THEN
-          MSQ=16*(AS+BS)*(CS+DS)*D03*D12/(S12-AMSS3(J)**2)**2
-         ELSE IF (PTYPE1.EQ.3.AND.KIN(3)) THEN
-          MSQ=16*(AS+BS)*(CS+DS)*D02*D13/(S13-AMSS3(J)**2)**2
-         ELSE IF (PTYPE1.EQ.4.AND.KIN(4)) THEN
-          MSQ=16*((AS+BS)*(CS+DS)*D01*D23+(AS-BS)*(CS+DS)*D23*
-     $        SQRT(AM0SQ*AM1SQ))/(S23-AMSS3(J)**2)**2
-         END IF
-        END IF          
-        IF (J.NE.JP) THEN
-        PTYPE2=KSS3(JP)
-        AP=ZISS3(1,JP)
-        BP=ZISS3(2,JP)
-        CP=ZFSS3(1,JP)
-        DP=ZFSS3(2,JP)
-        APC=CONJG(AP)
-        BPC=CONJG(BP)
-        CPC=CONJG(CP)
-        DPC=CONJG(DP)
-         IF (PTYPE1.EQ.2.AND.PTYPE2.EQ.2.AND.KIN(2).AND.KINP(2)) THEN
-          MMPD=16*D12*D03*(A*APC+B*BPC)*(C*CPC+D*DPC)/
-     $        (S12-AMSS3(J)**2)/(S12-AMSS3(JP)**2)
-          MSQ=2*REAL(MMPD)
-         END IF
-         IF (PTYPE1.EQ.3.AND.PTYPE2.EQ.3.AND.KIN(3).AND.KINP(3)) THEN
-          MMPD=16*D13*D02*(A*APC+B*BPC)*(C*CPC+D*DPC)/
-     $        (S13-AMSS3(J)**2)/(S13-AMSS3(JP)**2)
-          MSQ=2*REAL(MMPD)
-         END IF
-         IF (PTYPE1.EQ.4.AND.PTYPE2.EQ.4.AND.KIN(4).AND.KINP(4)) THEN
-          MMPD=16*D23*(D01*(A*APC+B*BPC)*(C*CPC+D*DPC)+
-     $        SQRT(AM0SQ*AM1SQ)*(A*APC-B*BPC)*(C*CPC-D*DPC))/
-     $        (S23-AMSS3(J)**2)/(S23-AMSS3(JP)**2)
-          MSQ=2*REAL(MMPD)
-         END IF
-         IF (PTYPE1.EQ.1.AND.PTYPE2.EQ.3.AND.KIN(1).AND.KINP(3)) THEN
-          MMPD=(16*D13*D02*((A*C-B*D)*(-APC*CPC+BPC*DPC)+
-     $         (A*D-B*C)*(APC*DPC-BPC*CPC))+
-     $         8*D23*SQRT(AM0SQ*AM1SQ)*((A*C+B*D)*(APC*CPC-BPC*DPC)-
-     $         (A*D+B*C)*(APC*DPC-BPC*CPC)))/
-     $         (S23-AMSS3(J)**2)/(S13-AMSS3(JP)**2)
-          MSQ=2*REAL(MMPD)
-         END IF
-         IF (PTYPE1.EQ.1.AND.PTYPE2.EQ.2.AND.KIN(1).AND.KINP(2)) THEN
-          MMPD=(16*D12*D03*((A*C+B*D)*(-APC*CPC+BPC*DPC)+
-     $         (A*D+B*C)*(APC*DPC-BPC*CPC))+
-     $         8*D23*SQRT(AM0SQ*AM1SQ)*((A*C-B*D)*(APC*CPC-BPC*DPC)+
-     $         (-A*D+B*C)*(APC*DPC+BPC*CPC)))/
-     $         (S23-AMSS3(J)**2)/(S12-AMSS3(JP)**2)
-          MSQ=2*REAL(MMPD)
-         END IF
-         IF (PTYPE1.EQ.3.AND.PTYPE2.EQ.4.AND.KIN(3).AND.KINP(4)) THEN
-          MMPD=((8*D13*D23+4*D23*AM1SQ)*((A*C+B*D)*(APC*CPC+BPC*DPC)+
-     $         (A*D+B*C)*(APC*DPC+BPC*CPC))+
-     $         4*D23*SQRT(AM0SQ*AM1SQ)*((A*C+B*D)*(APC*CPC-BPC*DPC)+
-     $         (A*D+B*C)*(APC*DPC-BPC*CPC)))/
-     $         (S13-AMSS3(J)**2)/(S23-AMSS3(JP)**2)
-          MSQ=2*REAL(MMPD)
-         END IF
-         IF (PTYPE1.EQ.2.AND.PTYPE2.EQ.4.AND.KIN(2).AND.KINP(4)) THEN
-          MMPD=-((8*D12*D23+4*D23*AM1SQ)*((A*C+B*D)*(APC*CPC+BPC*DPC)+
-     $         (A*D+B*C)*(APC*DPC+BPC*CPC))+
-     $         4*D23*SQRT(AM0SQ*AM1SQ)*((A*C+B*D)*(APC*CPC-BPC*DPC)+
-     $         (A*D+B*C)*(APC*DPC-BPC*CPC)))/
-     $         (S12-AMSS3(J)**2)/(S23-AMSS3(JP)**2)
-          MSQ=2*REAL(MMPD)
-         END IF
-         IF (PTYPE1.EQ.2.AND.PTYPE2.EQ.3.AND.KIN(2).AND.KINP(3)) THEN
-          MMPD=((8*D12*D13-4*D23*AM1SQ)*((A*C+B*D)*(APC*CPC+BPC*DPC)+
-     $         (A*D+B*C)*(APC*DPC+BPC*CPC))-
-     $         4*D23*SQRT(AM0SQ*AM1SQ)*((A*C-B*D)*(APC*CPC-BPC*DPC)+
-     $         (A*D-B*C)*(APC*DPC-BPC*CPC)))/
-     $         (S12-AMSS3(J)**2)/(S13-AMSS3(JP)**2)
-          MSQ=2*REAL(MMPD)
-         END IF
-        END IF
-        DECSS3=DECSS3+MSQ
-       END DO
-      END DO
-C
-      RETURN
-      END
diff --git a/ISAJET/code/dectau.F b/ISAJET/code/dectau.F
deleted file mode 100644 (file)
index b2b7325..0000000
+++ /dev/null
@@ -1,190 +0,0 @@
-#include "isajet/pilot.h"
-      LOGICAL FUNCTION DECTAU(IP,NADD,MEIP,IDABS,PREST)
-C
-C          Compute matrix elements for polarized tau decay.
-C          Polarization determined by tau parent.
-C          Auxiliary routine for DECAY. 
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-#include "isajet/itapes.inc"
-#include "isajet/wcon.inc"
-#include "isajet/partcl.inc"
-#include "isajet/dkytab.inc"
-#include "isajet/const.inc"
-#include "isajet/pjets.inc"
-#include "isajet/keys.inc"
-#include "isajet/xmssm.inc"
-#include "isajet/sspols.inc"
-#include "isajet/primar.inc"
-C
-      REAL PREST(4,6),WT,TAUHEL,S12,S12MAX,PIP,CTHNU,PSUM(4),AMV2,WT1
-      REAL DOT,DOT3,RANF,Z
-      INTEGER IP,NADD,IDABS(5),IPAR,IDPAR,JET,INU,I,K,I1,I2,IDSIB
-      INTEGER IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX,IDIP
-      INTEGER MEIP,IPX,IP1,IP2
-C
-      DOT(I1,I2)=PREST(4,I1)*PREST(4,I2)-PREST(1,I1)*PREST(1,I2)
-     $-PREST(2,I1)*PREST(2,I2)-PREST(3,I1)*PREST(3,I2)
-      DOT3(I1,I2)=PREST(1,I1)*PREST(1,I2)+PREST(2,I1)*PREST(2,I2)
-     $+PREST(3,I1)*PREST(3,I2)
-C
-      IDIP=IDENT(IP)
-      DECTAU=.TRUE.
-      IF(IABS(IDIP).NE.16) GO TO 999
-C
-C          Use PREST(K,6) for spin vector
-C
-      PIP=SQRT(PPTCL(1,IP)**2+PPTCL(2,IP)**2+PPTCL(3,IP)**2)
-      DO 100 K=1,3
-        PREST(K,6)=PPTCL(K,IP)/PIP
-100   CONTINUE
-      PREST(4,6)=0.
-C
-C          Take helicity TAUHEL=0 unless TAU parent is TP, W+-, H+-,
-C          or some SUSY particles.
-C          Take account of 1-particle decays!
-C
-
-      IPX=IP
-      TAUHEL=0.
-      IPAR=0
-      IDPAR=0
-110   IF(IORIG(IPX).GT.0) THEN
-        IPAR=MOD(IORIG(IPX),IPACK)
-        IDPAR=IDENT(IPAR)
-        IF(IDPAR.EQ.IDIP) THEN
-          IP1=IDCAY(IPAR)/IPACK
-          IP2=MOD(IDCAY(IPAR),IPACK)
-          IF(IP1.EQ.IP2) THEN
-            IPX=IPAR
-            GO TO 110
-          ENDIF
-        ENDIF
-        IDPAR=IABS(IDPAR)
-        IDSIB=0
-C          W/top parent
-        IF(IDPAR.GT.100.AND.MOD(IDPAR/10,10).GE.6) THEN
-          TAUHEL=-1.
-        ELSEIF(IDPAR.EQ.80) THEN
-          TAUHEL=-1.
-C          Charged Higgs parent
-        ELSEIF(IDPAR.EQ.86) THEN
-          TAUHEL=+1.
-C          SUSY parent - polarization also depends on sibling IDSIB
-        ELSEIF(GOMSSM.AND.IDPAR.GT.20.AND.IDPAR.LT.80) THEN
-          I1=IDCAY(IPAR)/IPACK
-          I2=MOD(IDCAY(IPAR),IPACK)
-          DO 120 I=I1,I2
-            IF(IABS(IDENT(I)).GT.20.AND.IABS(IDENT(I)).LT.80)
-     $      IDSIB=IABS(IDENT(I))
-120       CONTINUE
-          IF (IDPAR.EQ.35) THEN
-            TAUHEL=-1.
-          ELSEIF (IDPAR.EQ.36) THEN
-            IF (IDSIB.EQ.30) TAUHEL=PTAU1(1)
-            IF (IDSIB.EQ.40) TAUHEL=PTAU1(2)
-            IF (IDSIB.EQ.50) TAUHEL=PTAU1(3)
-            IF (IDSIB.EQ.60) TAUHEL=PTAU1(4)
-          ELSEIF (IDPAR.EQ.56) THEN
-            IF (IDSIB.EQ.30) TAUHEL=PTAU2(1)
-            IF (IDSIB.EQ.40) TAUHEL=PTAU2(2)
-            IF (IDSIB.EQ.50) TAUHEL=PTAU2(3)
-            IF (IDSIB.EQ.60) TAUHEL=PTAU2(4)
-          ELSEIF (IDPAR.EQ.39) THEN
-            IF(IDSIB.EQ.35) TAUHEL=-1.
-            IF(IDSIB.EQ.30) TAUHEL=PTAUWZ
-          ELSEIF (IDPAR.EQ.49.AND.IDSIB.EQ.35) THEN
-            TAUHEL=-1.
-          ELSEIF (IDPAR.EQ.40) THEN
-            IF(IDSIB.EQ.36) TAUHEL=PTAUZ2(1)
-            IF(IDSIB.EQ.56) TAUHEL=PTAUZ2(2)
-            IF(IDSIB.EQ.30) TAUHEL=PTAUZZ
-          ELSEIF (IDPAR.EQ.50) THEN
-            IF(IDSIB.EQ.36) TAUHEL=PTAUZ3(1)
-            IF(IDSIB.EQ.56) TAUHEL=PTAUZ3(2)
-          ELSEIF (IDPAR.EQ.60) THEN 
-            IF(IDSIB.EQ.36) TAUHEL=PTAUZ4(1)
-            IF(IDSIB.EQ.56) TAUHEL=PTAUZ4(2)
-          ENDIF
-        END IF
-      ELSE
-        IF(KEYS(3)) THEN
-          IF(IABS(IDENTW).EQ.80) TAUHEL=-1.
-        ELSE
-          JET=IABS(IORIG(IP))/IPACK
-          IF(JET.GT.0.AND.JET.LE.NJET) THEN
-            IF(IDJETS(JET).EQ.80) TAUHEL=-1.
-          ENDIF
-        ENDIF
-      ENDIF
-C
-C          Leptonic decays. DECTAU is always called for TAU- decay
-C          products, so selection is independent of IDENT(IP).
-C
-      IF(MEIP.EQ.5) THEN
-        IF(IDENT(NPTCL+1).LT.0) THEN
-          WT=PPTCL(5,IP)*(PREST(4,1)-TAUHEL*DOT(1,6))*DOT(2,3)
-        ELSEIF(IDENT(NPTCL+2).LT.0) THEN
-          WT=PPTCL(5,IP)*(PREST(4,2)-TAUHEL*DOT(2,6))*DOT(1,3)
-        ELSE
-          WT=PPTCL(5,IP)*(PREST(4,3)-TAUHEL*DOT(3,6))*DOT(1,2)
-        ENDIF
-        IF(WT.LT.RANF()*PPTCL(5,IP)**4/8.) THEN
-          DECTAU=.FALSE.
-        ELSE
-          DECTAU=.TRUE.
-        ENDIF
-        RETURN
-C
-C          Decay to PI + NUT, K + NUT
-C
-      ELSEIF(MEIP.EQ.6) THEN
-        INU=1
-        IF(IDABS(2).EQ.15) INU=2
-        CTHNU=DOT3(INU,6)/SQRT(DOT3(INU,INU))
-        WT=1.-TAUHEL*CTHNU
-        IF(WT.LT.RANF()*2.) THEN
-          DECTAU=.FALSE.
-        ELSE
-          DECTAU=.TRUE.
-        ENDIF
-        RETURN
-C
-C          Decay to RHO + NUT, A1 + NUT, K* + NUT
-C
-      ELSEIF(MEIP.EQ.7) THEN
-        DO 210 I=1,NADD
-210     IF(IDABS(I).EQ.15) INU=I
-        DO 220 K=1,4
-          PSUM(K)=0.
-          DO 221 I=1,NADD
-            IF(I.EQ.INU) GO TO 221
-            PSUM(K)=PSUM(K)+PREST(K,I)
-221       CONTINUE
-220     CONTINUE
-        AMV2=PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2
-        WT1=2.*AMV2/(2.*AMV2+PPTCL(5,IP)**2)
-        CTHNU=DOT3(INU,6)/SQRT(DOT3(INU,INU))
-        WT=WT1*(1.+TAUHEL*CTHNU)+(1.-WT1)*(1-TAUHEL*CTHNU)
-        IF(WT.LT.RANF()*2.) THEN
-          DECTAU=.FALSE.
-        ELSE
-          DECTAU=.TRUE.
-        ENDIF
-        RETURN
-C
-C          Ignore matrix element for all other decays
-C
-      ELSE
-        DECTAU=.TRUE.
-        RETURN
-      ENDIF
-C
-C          Error
-C
-999   CALL PRTEVT(0)
-      WRITE(ITLIS,99999) IP
-99999 FORMAT(//5X,'ERROR IN DECTAU FOR PARTICLE',I6)
-      END
diff --git a/ISAJET/code/decva.F b/ISAJET/code/decva.F
deleted file mode 100644 (file)
index 64b9cf3..0000000
+++ /dev/null
@@ -1,44 +0,0 @@
-#include "isajet/pilot.h"
-      LOGICAL FUNCTION DECVA(IP,NADD,IDABS,PREST)
-C
-C          Compute matrix element unpolarized for V-A.
-C          Auxiliary routine for DECAY. 
-C
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-#include "isajet/itapes.inc"
-#include "isajet/wcon.inc"
-#include "isajet/partcl.inc"
-#include "isajet/dkytab.inc"
-#include "isajet/jetset.inc"
-#include "isajet/jwork.inc"
-#include "isajet/const.inc"
-#include "isajet/keys.inc"
-#include "isajet/pjets.inc"
-#include "isajet/xmssm.inc"
-#include "isajet/sspols.inc"
-C
-      REAL PREST(4,6)
-      REAL DOT,RANF,WT
-      INTEGER IP,NADD,IDABS(5),I,K,I1,I2,IDIPA
-C
-      DOT(I1,I2)=PREST(4,I1)*PREST(4,I2)-PREST(1,I1)*PREST(1,I2)
-     $-PREST(2,I1)*PREST(2,I2)-PREST(3,I1)*PREST(3,I2)
-C
-      IDIPA=IABS(IDENT(IP))
-      IF(IDENT(NPTCL+1).LT.0) THEN
-        WT=PPTCL(5,IP)*PREST(4,1)*DOT(2,3)
-      ELSEIF(IDENT(NPTCL+2).LT.0) THEN
-        WT=PPTCL(5,IP)*PREST(4,2)*DOT(1,3)
-      ELSE
-        WT=PPTCL(5,IP)*PREST(4,3)*DOT(1,2)
-      ENDIF
-      IF(WT.LT.RANF()*PPTCL(5,IP)**4/16.) THEN
-        DECVA=.FALSE.
-      ELSE
-        DECVA=.TRUE.
-      ENDIF
-      RETURN
-      END
diff --git a/ISAJET/code/dhelas.F b/ISAJET/code/dhelas.F
deleted file mode 100644 (file)
index 17058af..0000000
+++ /dev/null
@@ -1,3748 +0,0 @@
-#include "isajet/pilot.h"
-C  *********************************************************************
-C  ***                                                               ***
-C  ***               coded by H. Murayama & I. Watanabe              ***
-C  *** For the formalism and notations, see the following reference: ***
-C  ***           H. Murayama, I. Watanabe and K. Hagiwara            ***
-C  ***           "HELAS: HELicity Amplitude Subroutines              ***
-C  ***               for Feynman diagram evaluation"                 ***
-C  ***               KEK Report 91-11, December 1991                 ***
-C  ***                                                               ***
-C  *********************************************************************
-C
-C  Converted to double precision by W. Long and T. Seltzer for MadGraph.
-C
-C  Minor changes for portability by FEP, July 1999. The code is not ANSI
-C  standard, but that cannot be helped if MadGraph compatibility is to 
-C  be maintained.
-C
-C ======================================================================
-C
-      SUBROUTINE BOOSTX(P,Q , PBOOST)
-C
-C this subroutine performs the lorentz boost of a four-momentum.  the   
-C momentum p is assumed to be given in the rest frame of q.  pboost is  
-C the momentum p boosted to the frame in which q is given.  q must be a 
-C timelike momentum.                                                    
-C                                                                       
-C input:                                                                
-C       real    p(0:3)         : four-momentum p in the q rest  frame   
-C       real    q(0:3)         : four-momentum q in the boosted frame   
-C                                                                       
-C output:                                                               
-C       real    pboost(0:3)    : four-momentum p in the boosted frame   
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      REAL*8    P(0:3),Q(0:3),PBOOST(0:3),PQ,QQ,M,LF
-      REAL*8 RXZERO
-      PARAMETER( RXZERO=0.0D0 )
-C
-      QQ=Q(1)**2+Q(2)**2+Q(3)**2
-C
-      IF ( QQ .NE. RXZERO ) THEN
-         PQ=P(1)*Q(1)+P(2)*Q(2)+P(3)*Q(3)
-         M=SQRT(Q(0)**2-QQ)
-         LF=((Q(0)-M)*PQ/QQ+P(0))/M
-         PBOOST(0) = (P(0)*Q(0)+PQ)/M
-         PBOOST(1) =  P(1)+Q(1)*LF
-         PBOOST(2) =  P(2)+Q(2)*LF
-         PBOOST(3) =  P(3)+Q(3)*LF
-      ELSE
-         PBOOST(0)=P(0)
-         PBOOST(1)=P(1)
-         PBOOST(2)=P(2)
-         PBOOST(3)=P(3)
-      ENDIF
-C
-      RETURN
-      END
-C
-C **********************************************************************
-C
-      SUBROUTINE COUP1X(SW2 , GW,GWWA,GWWZ)
-C
-C this subroutine sets up the coupling constants of the gauge bosons in 
-C the standard model.                                                   
-C                                                                       
-C input:                                                                
-C       real    sw2            : square of sine of the weak angle       
-C                                                                       
-C output:                                                               
-C       real    gw             : weak coupling constant                 
-C       real    gwwa           : dimensionless coupling of w-,w+,a      
-C       real    gwwz           : dimensionless coupling of w-,w+,z      
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      REAL*8    SW2,GW,GWWA,GWWZ,ALPHA,FOURPI,EE,SW,CW
-      REAL*8 RXONE, RXFOUR, RXOTE, RXPI, RIALPH
-      PARAMETER( RXONE=1.0D0, RXFOUR=4.0D0, RXOTE=128.0D0 )
-      PARAMETER( RXPI=3.14159265358979323846D0, RIALPH=137.0359895D0 )
-C
-      ALPHA = RXONE / RXOTE
-C      alpha = r_one / r_ialph
-      FOURPI = RXFOUR * RXPI
-      EE=SQRT( ALPHA * FOURPI )
-      SW=SQRT( SW2 )
-      CW=SQRT( RXONE - SW2 )
-C
-      GW    =  EE/SW
-      GWWA  =  EE
-      GWWZ  =  EE*CW/SW
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE COUP2X(SW2 , GAL,GAU,GAD,GWF,GZN,GZL,GZU,GZD,G1)
-C
-C this subroutine sets up the coupling constants for the fermion-       
-C fermion-vector vertices in the standard model.  the array of the      
-C couplings specifies the chirality of the flowing-in fermion.  g??(1)  
-C denotes a left-handed coupling, and g??(2) a right-handed coupling.   
-C                                                                       
-C input:                                                                
-C       real    sw2            : square of sine of the weak angle       
-C                                                                       
-C output:                                                               
-C       real    gal(2)         : coupling with a of charged leptons     
-C       real    gau(2)         : coupling with a of up-type quarks      
-C       real    gad(2)         : coupling with a of down-type quarks    
-C       real    gwf(2)         : coupling with w-,w+ of fermions        
-C       real    gzn(2)         : coupling with z of neutrinos           
-C       real    gzl(2)         : coupling with z of charged leptons     
-C       real    gzu(2)         : coupling with z of up-type quarks      
-C       real    gzd(2)         : coupling with z of down-type quarks    
-C       real    g1(2)          : unit coupling of fermions              
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      REAL*8 GAL(2),GAU(2),GAD(2),GWF(2),GZN(2),GZL(2),GZU(2),GZD(2),
-     &     G1(2),SW2,ALPHA,FOURPI,EE,SW,CW,EZ,EY
-C
-      REAL*8 RXZERO, RXHALF, RXONE, RXTWO, RTHREE, RXFOUR, RXOTE
-      REAL*8 RXPI, RIALPH
-      PARAMETER( RXZERO=0.0D0, RXHALF=0.5D0, RXONE=1.0D0, RXTWO=2.0D0,
-     $     RTHREE=3.0D0 )
-      PARAMETER( RXFOUR=4.0D0, RXOTE=128.0D0 )
-      PARAMETER( RXPI=3.14159265358979323846D0, RIALPH=137.0359895D0 )
-C
-      ALPHA = RXONE / RXOTE
-C      alpha = r_one / r_ialph
-      FOURPI = RXFOUR * RXPI
-      EE=SQRT( ALPHA * FOURPI )
-      SW=SQRT( SW2 )
-      CW=SQRT( RXONE - SW2 )
-      EZ=EE/(SW*CW)
-      EY=EE*(SW/CW)
-C
-      GAL(1) =  EE
-      GAL(2) =  EE
-      GAU(1) = -EE*RXTWO/RTHREE
-      GAU(2) = -EE*RXTWO/RTHREE
-      GAD(1) =  EE   /RTHREE
-      GAD(2) =  EE   /RTHREE
-      GWF(1) = -EE/SQRT(RXTWO*SW2)
-      GWF(2) =  RXZERO
-      GZN(1) = -EZ*  RXHALF
-      GZN(2) =  RXZERO
-      GZL(1) = -EZ*(-RXHALF+SW2)
-      GZL(2) = -EY
-      GZU(1) = -EZ*( RXHALF-SW2*RXTWO/RTHREE)
-      GZU(2) =  EY*          RXTWO/RTHREE
-      GZD(1) = -EZ*(-RXHALF+SW2   /RTHREE)
-      GZD(2) = -EY             /RTHREE
-      G1(1)  =  RXONE
-      G1(2)  =  RXONE
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE COUP3X(SW2,ZMASS,HMASS , 
-     &                  GWWH,GZZH,GHHH,GWWHH,GZZHH,GHHHH)
-C
-C this subroutine sets up the coupling constants of the gauge bosons and
-C higgs boson in the standard model.                                    
-C                                                                       
-C input:                                                                
-C       real    sw2            : square of sine of the weak angle       
-C       real    zmass          : mass of z                              
-C       real    hmass          : mass of higgs                          
-C                                                                       
-C output:                                                               
-C       real    gwwh           : dimensionful  coupling of w-,w+,h      
-C       real    gzzh           : dimensionful  coupling of z, z, h      
-C       real    ghhh           : dimensionful  coupling of h, h, h      
-C       real    gwwhh          : dimensionful  coupling of w-,w+,h, h   
-C       real    gzzhh          : dimensionful  coupling of z, z, h, h   
-C       real    ghhhh          : dimensionless coupling of h, h, h, h   
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      REAL*8    SW2,ZMASS,HMASS,GWWH,GZZH,GHHH,GWWHH,GZZHH,GHHHH,
-     &        ALPHA,FOURPI,EE2,SC2,V
-C
-      REAL*8 RXHALF, RXONE, RXTWO, RTHREE, RXFOUR, RXOTE
-      REAL*8 RXPI, RIALPH
-      PARAMETER( RXHALF=0.5D0, RXONE=1.0D0, RXTWO=2.0D0, RTHREE=3.0D0 )
-      PARAMETER( RXFOUR=4.0D0, RXOTE=128.0D0 )
-      PARAMETER( RXPI=3.14159265358979323846D0, RIALPH=137.0359895D0 )
-C
-      ALPHA = RXONE / RXOTE
-C      alpha = r_one / r_ialph
-      FOURPI = RXFOUR * RXPI
-      EE2=ALPHA*FOURPI
-      SC2=SW2*( RXONE - SW2 )
-      V = RXTWO * ZMASS*SQRT(SC2)/SQRT(EE2)
-C
-      GWWH  =   EE2/SW2*RXHALF*V
-      GZZH  =   EE2/SC2*RXHALF*V
-      GHHH  =  -HMASS**2/V*RTHREE
-      GWWHH =   EE2/SW2*RXHALF
-      GZZHH =   EE2/SC2*RXHALF
-      GHHHH = -(HMASS/V)**2*RTHREE
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE COUP4X(SW2,ZMASS,FMASS , GCHF)
-C
-C This subroutine sets up the coupling constant for the fermion-fermion-
-C Higgs vertex in the STANDARD MODEL.  The coupling is COMPLEX and the  
-C array of the coupling specifies the chirality of the flowing-IN       
-C fermion.  GCHF(1) denotes a left-handed coupling, and GCHF(2) a right-
-C handed coupling.                                                      
-C                                                                       
-C INPUT:                                                                
-C       real    SW2            : square of sine of the weak angle       
-C       real    ZMASS          : Z       mass                           
-C       real    FMASS          : fermion mass                           
-C                                                                       
-C OUTPUT:                                                               
-C       complex GCHF(2)        : coupling of fermion and Higgs          
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 GCHF(2)
-      REAL*8    SW2,ZMASS,FMASS,ALPHA,FOURPI,EZ,G
-C
-      ALPHA=1.D0/128.D0
-C      ALPHA=1./REAL(137.0359895)
-      FOURPI=4.D0*3.14159265358979323846D0
-      EZ=SQRT(ALPHA*FOURPI)/SQRT(SW2*(1.D0-SW2))
-      G=EZ*FMASS*0.5D0/ZMASS
-C
-      GCHF(1) = DCMPLX( -G )
-      GCHF(2) = DCMPLX( -G )
-C
-      RETURN
-      END
-C
-C ======================================================================
-C
-      SUBROUTINE EAIXXX(EB,EA,SHLF,CHLF,PHI,NHE,NHA , EAI)
-C
-C This subroutine computes an off-shell electron wavefunction after     
-C emitting a photon from the electron beam, with a special care for the 
-C small angle region.  The momenta are measured in the laboratory frame,
-C where the e- beam is along the positive z axis.                       
-C                                                                       
-C INPUT:                                                                
-C       real    EB             : energy (GeV)    of beam  e-            
-C       real    EA             : energy (GeV)    of final photon        
-C       real    SHLF           : sin(theta/2)    of final photon        
-C       real    CHLF           : cos(theta/2)    of final photon        
-C       real    PHI            : azimuthal angle of final photon        
-C       integer NHE  = -1 or 1 : helicity        of beam  e-            
-C       integer NHA  = -1 or 1 : helicity        of final photon        
-C                                                                       
-C OUTPUT:                                                               
-C       complex EAI(6)         : off-shell electron             |e',A,e>
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 EAI(6),PHS
-      REAL*8  EB,EA,SHLF,CHLF,PHI,ME,ALPHA,GAL,RNHE,X,C,S,D,COEFF,
-     &        XNNP,XNNM,SNP,CSP
-      INTEGER NHE,NHA,NN
-C
-      ME   = 0.51099906D-3
-      ALPHA=1./128.
-      GAL  =SQRT(ALPHA*4.*3.14159265D0)
-C
-      NN=NHA*NHE
-      RNHE=NHE
-      X=EA/EB
-      C=(CHLF+SHLF)*(CHLF-SHLF)
-      S=2.*CHLF*SHLF
-      D=-1./(EA*EB*(4.*SHLF**2+(ME/EB)**2*C))
-      COEFF=-NN*GAL*SQRT(EB)*D
-      XNNP=X*(1+NN)
-      XNNM=X*(1-NN)
-      SNP=SIN(PHI)
-      CSP=COS(PHI)
-      PHS=DCMPLX( CSP , RNHE*SNP )
-C
-      EAI((5-3*NHE)/2) = -RNHE*COEFF*ME*S*(1.+XNNP*.5)
-      EAI((5-NHE)/2)   =  XNNP*COEFF*ME*CHLF**2*PHS
-      EAI((5+NHE)/2)   =  RNHE*COEFF*EB*S*(-2.+XNNM)
-      EAI((5+3*NHE)/2) =  XNNM*COEFF*EB*SHLF**2*PHS*2.
-C
-      EAI(5) =  EB*DCMPLX( 1.-X , 1.-X*C )
-      EAI(6) = -EB*X*S*DCMPLX( CSP , SNP )
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE EAOXXX(EB,EA,SHLF,CHLF,PHI,NHE,NHA , EAO)
-C
-C This subroutine computes an off-shell positron wavefunction after     
-C emitting a photon from the positron beam, with a special care for the 
-C small angle region.  The momenta are measured in the laboratory frame,
-C where the e+ beam is along the negative z axis.                       
-C                                                                       
-C INPUT:                                                                
-C       real    EB             : energy (GeV)    of beam  e+            
-C       real    EA             : energy (GeV)    of final photon        
-C       real    SHLF           : sin(theta/2)    of final photon        
-C       real    CHLF           : cos(theta/2)    of final photon        
-C       real    PHI            : azimuthal angle of final photon        
-C       integer NHE  = -1 or 1 : helicity        of beam  e+            
-C       integer NHA  = -1 or 1 : helicity        of final photon        
-C                                                                       
-C OUTPUT:                                                               
-C       complex EAO(6)         : off-shell positron             <e,A,e'|
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 EAO(6),PHS
-      REAL*8  EB,EA,SHLF,CHLF,PHI,ME,ALPHA,GAL,RNHE,X,C,S,D,COEFF,
-     &        XNNP,XNNM,SNP,CSP
-      INTEGER NHE,NHA,NN
-C
-      ME   = 0.51099906D-3
-      ALPHA=1./128.
-      GAL  =SQRT(ALPHA*4.*3.14159265D0)
-C
-      NN=NHA*NHE
-      RNHE=NHE
-      X=EA/EB
-      C=(CHLF+SHLF)*(CHLF-SHLF)
-      S=2.*CHLF*SHLF
-      D=-1./(EA*EB*(4.*CHLF**2-(ME/EB)**2*C))
-      COEFF=NN*GAL*SQRT(EB)*D
-      XNNP=X*(1+NN)
-      XNNM=X*(1-NN)
-      SNP=SIN(PHI)
-      CSP=COS(PHI)
-      PHS=DCMPLX( CSP ,-RNHE*SNP )
-C
-      EAO((5-3*NHE)/2) =              COEFF*ME*S*(1.+XNNP*.5)
-      EAO((5-NHE)/2)   = RNHE*XNNP    *COEFF*ME*SHLF**2*PHS
-      EAO((5+NHE)/2)   =              COEFF*EB*S*(-2.+XNNM)
-      EAO((5+3*NHE)/2) = REAL(NHA-NHE)*COEFF*EB*X*CHLF**2*PHS*2.
-C
-      EAO(5) = EB*DCMPLX( X-1. , X*C+1. )
-      EAO(6) = EB*X*S*DCMPLX( CSP , SNP )
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE FSIXXX(FI,SC,GC,FMASS,FWIDTH , FSI)
-C
-C this subroutine computes an off-shell fermion wavefunction from a     
-C flowing-in external fermion and a vector boson.                       
-C                                                                       
-C input:                                                                
-C       complex*16 fi(6)          : flow-in  fermion                   |fi>
-C       complex*16 sc(3)          : input    scalar                      s 
-C       complex*16 gc(2)          : coupling constants                 gchf
-C       real*8    fmass          : mass  of output fermion f'             
-C       real*8    fwidth         : width of output fermion f'             
-C                                                                       
-C output:                                                               
-C       complex fsi(6)         : off-shell fermion             |f',s,fi>
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FI(6),SC(3),FSI(6),GC(2),SL1,SL2,SR1,SR2,DS
-      REAL*8     PF(0:3),FMASS,FWIDTH,PF2,P0P3,P0M3
-C
-      FSI(5) = FI(5)-SC(2)
-      FSI(6) = FI(6)-SC(3)
-C
-      PF(0)=DBLE( FSI(5))
-      PF(1)=DBLE( FSI(6))
-      PF(2)=DIMAG(FSI(6))
-      PF(3)=DIMAG(FSI(5))
-      PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
-C
-      DS=-SC(1)/DCMPLX(PF2-FMASS**2,MAX(DSIGN(FMASS*FWIDTH ,PF2),0D0))
-      P0P3=PF(0)+PF(3)
-      P0M3=PF(0)-PF(3)
-      SL1=GC(1)*(P0P3*FI(1)+DCONJG(FSI(6))*FI(2))
-      SL2=GC(1)*(P0M3*FI(2)      +FSI(6) *FI(1))
-      SR1=GC(2)*(P0M3*FI(3)-DCONJG(FSI(6))*FI(4))
-      SR2=GC(2)*(P0P3*FI(4)      -FSI(6) *FI(3))
-C
-      FSI(1) = ( GC(1)*FMASS*FI(1) + SR1 )*DS
-      FSI(2) = ( GC(1)*FMASS*FI(2) + SR2 )*DS
-      FSI(3) = ( GC(2)*FMASS*FI(3) + SL1 )*DS
-      FSI(4) = ( GC(2)*FMASS*FI(4) + SL2 )*DS
-C
-      RETURN          
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE FSOXXX(FO,SC,GC,FMASS,FWIDTH , FSO)
-C
-C this subroutine computes an off-shell fermion wavefunction from a     
-C flowing-out external fermion and a vector boson.                      
-C                                                                       
-C input:                                                                
-C       complex*16 fo(6)          : flow-out fermion                   <fo|
-C       complex*16 sc(6)          : input    scalar                      s 
-C       complex*16 gc(2)          : coupling constants                 gchf
-C       real*8     fmass          : mass  of output fermion f'             
-C       real*8     fwidth         : width of output fermion f'             
-C                                                                       
-C output:                                                               
-C       complex fso(6)         : off-shell fermion             <fo,s,f'|
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FO(6),SC(6),FSO(6),GC(2),SL1,SL2,SR1,SR2,DS
-      REAL*8     PF(0:3),FMASS,FWIDTH,PF2,P0P3,P0M3
-C
-      FSO(5) = FO(5)+SC(2)
-      FSO(6) = FO(6)+SC(3)
-C
-      PF(0)=DBLE( FSO(5))
-      PF(1)=DBLE( FSO(6))
-      PF(2)=DIMAG(FSO(6))
-      PF(3)=DIMAG(FSO(5))
-      PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
-C
-      DS=-SC(1)/DCMPLX(PF2-FMASS**2,MAX(DSIGN(FMASS*FWIDTH ,PF2),0D0))
-      P0P3=PF(0)+PF(3)
-      P0M3=PF(0)-PF(3)
-      SL1=GC(2)*(P0P3*FO(3)      +FSO(6) *FO(4))
-      SL2=GC(2)*(P0M3*FO(4)+DCONJG(FSO(6))*FO(3))
-      SR1=GC(1)*(P0M3*FO(1)      -FSO(6) *FO(2))
-      SR2=GC(1)*(P0P3*FO(2)-DCONJG(FSO(6))*FO(1))
-C
-      FSO(1) = ( GC(1)*FMASS*FO(1) + SL1 )*DS
-      FSO(2) = ( GC(1)*FMASS*FO(2) + SL2 )*DS
-      FSO(3) = ( GC(2)*FMASS*FO(3) + SR1 )*DS
-      FSO(4) = ( GC(2)*FMASS*FO(4) + SR2 )*DS
-C
-      RETURN          
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE FVIXXX(FI,VC,G,FMASS,FWIDTH , FVI)
-C
-C this subroutine computes an off-shell fermion wavefunction from a     
-C flowing-in external fermion and a vector boson.                       
-C                                                                       
-C input:                                                                
-C       complex fi(6)          : flow-in  fermion                   |fi>
-C       complex vc(6)          : input    vector                      v 
-C       real    g(2)           : coupling constants                  gvf
-C       real    fmass          : mass  of output fermion f'             
-C       real    fwidth         : width of output fermion f'             
-C                                                                       
-C output:                                                               
-C       complex fvi(6)         : off-shell fermion             |f',v,fi>
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FI(6),VC(6),FVI(6),SL1,SL2,SR1,SR2,D
-      REAL*8    G(2),PF(0:3),FMASS,FWIDTH,PF2
-      REAL*8 RXZERO, RXONE
-      PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
-      COMPLEX*16 CXIMAG
-C
-      LOGICAL FIRST
-      SAVE CXIMAG,FIRST
-      DATA FIRST/.TRUE./
-C
-C          Fix compilation with g77
-      IF(FIRST) THEN
-        FIRST=.FALSE.
-        CXIMAG=DCMPLX( RXZERO, RXONE )
-      ENDIF
-C
-      FVI(5) = FI(5)-VC(5)
-      FVI(6) = FI(6)-VC(6)
-C
-      PF(0)=DBLE( FVI(5))
-      PF(1)=DBLE( FVI(6))
-      PF(2)=DIMAG(FVI(6))
-      PF(3)=DIMAG(FVI(5))
-      PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
-C
-      D=-RXONE/DCMPLX( PF2-FMASS**2,MAX(SIGN(FMASS*FWIDTH,PF2),RXZERO))
-      SL1= (VC(1)+       VC(4))*FI(1)
-     &    +(VC(2)-CXIMAG*VC(3))*FI(2)
-      SL2= (VC(2)+CXIMAG*VC(3))*FI(1)
-     &    +(VC(1)-       VC(4))*FI(2)
-C
-      IF ( G(2) .NE. RXZERO ) THEN
-         SR1= (VC(1)-       VC(4))*FI(3)
-     &       -(VC(2)-CXIMAG*VC(3))*FI(4)
-         SR2=-(VC(2)+CXIMAG*VC(3))*FI(3)
-     &       +(VC(1)+       VC(4))*FI(4)
-C
-         FVI(1) = ( G(1)*((PF(0)-PF(3))*SL1 -DCONJG(FVI(6))*SL2)
-     &             +G(2)*FMASS*SR1)*D
-         FVI(2) = ( G(1)*(      -FVI(6)*SL1 +(PF(0)+PF(3))*SL2)
-     &             +G(2)*FMASS*SR2)*D
-         FVI(3) = ( G(2)*((PF(0)+PF(3))*SR1 +DCONJG(FVI(6))*SR2)
-     &             +G(1)*FMASS*SL1)*D
-         FVI(4) = ( G(2)*(       FVI(6)*SR1 +(PF(0)-PF(3))*SR2)
-     &             +G(1)*FMASS*SL2)*D
-C
-      ELSE          
-         FVI(1) = G(1)*((PF(0)-PF(3))*SL1 -DCONJG(FVI(6))*SL2)*D
-         FVI(2) = G(1)*(      -FVI(6)*SL1 +(PF(0)+PF(3))*SL2)*D
-         FVI(3) = G(1)*FMASS*SL1*D
-         FVI(4) = G(1)*FMASS*SL2*D
-      END IF
-C
-      RETURN          
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE FVOXXX(FO,VC,G,FMASS,FWIDTH , FVO)
-C
-C this subroutine computes an off-shell fermion wavefunction from a     
-C flowing-out external fermion and a vector boson.                      
-C                                                                       
-C input:                                                                
-C       complex fo(6)          : flow-out fermion                   <fo|
-C       complex vc(6)          : input    vector                      v 
-C       real    g(2)           : coupling constants                  gvf
-C       real    fmass          : mass  of output fermion f'             
-C       real    fwidth         : width of output fermion f'             
-C                                                                       
-C output:                                                               
-C       complex fvo(6)         : off-shell fermion             <fo,v,f'|
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FO(6),VC(6),FVO(6),SL1,SL2,SR1,SR2,D
-      REAL*8    G(2),PF(0:3),FMASS,FWIDTH,PF2
-      REAL*8 RXZERO, RXONE
-      PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
-      COMPLEX*16 CXIMAG
-      LOGICAL FIRST
-      SAVE CXIMAG,FIRST
-      DATA FIRST/.TRUE./
-C
-C          Fix compilation with g77
-      IF(FIRST) THEN
-        FIRST=.FALSE.
-        CXIMAG=DCMPLX( RXZERO, RXONE )
-      ENDIF
-C
-      FVO(5) = FO(5)+VC(5)
-      FVO(6) = FO(6)+VC(6)
-C
-      PF(0)=DBLE( FVO(5))
-      PF(1)=DBLE( FVO(6))
-      PF(2)=DIMAG(FVO(6))
-      PF(3)=DIMAG(FVO(5))
-      PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
-C
-      D=-RXONE/DCMPLX( PF2-FMASS**2,MAX(SIGN(FMASS*FWIDTH,PF2),RXZERO))
-      SL1= (VC(1)+       VC(4))*FO(3)
-     &    +(VC(2)+CXIMAG*VC(3))*FO(4)
-      SL2= (VC(2)-CXIMAG*VC(3))*FO(3)
-     &    +(VC(1)-       VC(4))*FO(4)
-C
-      IF ( G(2) .NE. RXZERO ) THEN
-         SR1= (VC(1)-       VC(4))*FO(1)
-     &       -(VC(2)+CXIMAG*VC(3))*FO(2)
-         SR2=-(VC(2)-CXIMAG*VC(3))*FO(1)
-     &       +(VC(1)+       VC(4))*FO(2)
-C
-         FVO(1) = ( G(2)*( (PF(0)+PF(3))*SR1        +FVO(6)*SR2)
-     &             +G(1)*FMASS*SL1)*D
-         FVO(2) = ( G(2)*( DCONJG(FVO(6))*SR1 +(PF(0)-PF(3))*SR2)
-     &             +G(1)*FMASS*SL2)*D
-         FVO(3) = ( G(1)*( (PF(0)-PF(3))*SL1        -FVO(6)*SL2)
-     &             +G(2)*FMASS*SR1)*D
-         FVO(4) = ( G(1)*(-DCONJG(FVO(6))*SL1 +(PF(0)+PF(3))*SL2)
-     &             +G(2)*FMASS*SR2)*D
-C
-      ELSE
-         FVO(1) = G(1)*FMASS*SL1*D
-         FVO(2) = G(1)*FMASS*SL2*D
-         FVO(3) = G(1)*( (PF(0)-PF(3))*SL1        -FVO(6)*SL2)*D
-         FVO(4) = G(1)*(-DCONJG(FVO(6))*SL1 +(PF(0)+PF(3))*SL2)*D
-      END IF
-C
-      RETURN          
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE GGGGXX(WM,W31,WP,W32,G, VERTEX)
-C
-C this subroutine computes an amplitude of the four-point coupling of   
-C the w-, w+ and two w3/z/a.  the amplitude includes the contributions  
-C of w exchange diagrams.  the internal w propagator is given in unitary
-C gauge.  if one sets wmass=0.0, then the gggg vertex is given (see sect
-C 2.9.1 of the manual).
-C                                                                       
-C input:                                                                
-C       complex wm(0:3)        : flow-out w-                         wm 
-C       complex w31(0:3)       : first    w3/z/a                     w31
-C       complex wp(0:3)        : flow-out w+                         wp 
-C       complex w32(0:3)       : second   w3/z/a                     w32
-C       real    g              : coupling of w31 with w-/w+             
-C                                                  (see the table below)
-C                                                                       
-C the possible sets of the inputs are as follows:                       
-C   -------------------------------------------                         
-C   |  wm  |  w31 |  wp  |  w32 |  g31 |  g32 |                         
-C   -------------------------------------------                         
-C   |  w-  |  w3  |  w+  |  w3  |  gw  |  gw  |                         
-C   |  w-  |  w3  |  w+  |  z   |  gw  | gwwz |                         
-C   |  w-  |  w3  |  w+  |  a   |  gw  | gwwa |                         
-C   |  w-  |  z   |  w+  |  z   | gwwz | gwwz |                         
-C   |  w-  |  z   |  w+  |  a   | gwwz | gwwa |                         
-C   |  w-  |  a   |  w+  |  a   | gwwa | gwwa |                         
-C   -------------------------------------------                         
-C where all the bosons are defined by the flowing-out quantum number.   
-C                                                                       
-C output:                                                               
-C       complex vertex         : amplitude          gamma(wm,w31,wp,w32)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16    WM(6),W31(6),WP(6),W32(6),VERTEX
-      COMPLEX*16 DV1(0:3),DV2(0:3),DV3(0:3),DV4(0:3),
-     &           DVERTX,V12,V13,V14,V23,V24,V34
-      REAL*8       PWM(0:3),PW31(0:3),PWP(0:3),PW32(0:3),G
-      REAL*8     DP1(0:3),DP2(0:3),DP3(0:3),DP4(0:3)
-      REAL*8 RXZERO, RXONE
-      PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
-C
-      PWM(0)=DBLE( WM(5))
-      PWM(1)=DBLE( WM(6))
-      PWM(2)=DIMAG(WM(6))
-      PWM(3)=DIMAG(WM(5))
-      PWP(0)=DBLE( WP(5))
-      PWP(1)=DBLE( WP(6))
-      PWP(2)=DIMAG(WP(6))
-      PWP(3)=DIMAG(WP(5))
-      PW31(0)=DBLE( W31(5))
-      PW31(1)=DBLE( W31(6))
-      PW31(2)=DIMAG(W31(6))
-      PW31(3)=DIMAG(W31(5))
-      PW32(0)=DBLE( W32(5))
-      PW32(1)=DBLE( W32(6))
-      PW32(2)=DIMAG(W32(6))
-      PW32(3)=DIMAG(W32(5))
-C
-      DV1(0)=DCMPLX(WM(1))
-      DV1(1)=DCMPLX(WM(2))
-      DV1(2)=DCMPLX(WM(3))
-      DV1(3)=DCMPLX(WM(4))
-      DP1(0)=DBLE(PWM(0))
-      DP1(1)=DBLE(PWM(1))
-      DP1(2)=DBLE(PWM(2))
-      DP1(3)=DBLE(PWM(3))
-      DV2(0)=DCMPLX(W31(1))
-      DV2(1)=DCMPLX(W31(2))
-      DV2(2)=DCMPLX(W31(3))
-      DV2(3)=DCMPLX(W31(4))
-      DP2(0)=DBLE(PW31(0))
-      DP2(1)=DBLE(PW31(1))
-      DP2(2)=DBLE(PW31(2))
-      DP2(3)=DBLE(PW31(3))
-      DV3(0)=DCMPLX(WP(1))
-      DV3(1)=DCMPLX(WP(2))
-      DV3(2)=DCMPLX(WP(3))
-      DV3(3)=DCMPLX(WP(4))
-      DP3(0)=DBLE(PWP(0))
-      DP3(1)=DBLE(PWP(1))
-      DP3(2)=DBLE(PWP(2))
-      DP3(3)=DBLE(PWP(3))
-      DV4(0)=DCMPLX(W32(1))
-      DV4(1)=DCMPLX(W32(2))
-      DV4(2)=DCMPLX(W32(3))
-      DV4(3)=DCMPLX(W32(4))
-      DP4(0)=DBLE(PW32(0))
-      DP4(1)=DBLE(PW32(1))
-      DP4(2)=DBLE(PW32(2))
-      DP4(3)=DBLE(PW32(3))
-C
-      V12= DV1(0)*DV2(0)-DV1(1)*DV2(1)-DV1(2)*DV2(2)-DV1(3)*DV2(3)
-      V13= DV1(0)*DV3(0)-DV1(1)*DV3(1)-DV1(2)*DV3(2)-DV1(3)*DV3(3)
-      V14= DV1(0)*DV4(0)-DV1(1)*DV4(1)-DV1(2)*DV4(2)-DV1(3)*DV4(3)
-      V23= DV2(0)*DV3(0)-DV2(1)*DV3(1)-DV2(2)*DV3(2)-DV2(3)*DV3(3)
-      V24= DV2(0)*DV4(0)-DV2(1)*DV4(1)-DV2(2)*DV4(2)-DV2(3)*DV4(3)
-      V34= DV3(0)*DV4(0)-DV3(1)*DV4(1)-DV3(2)*DV4(2)-DV3(3)*DV4(3)
-
-         DVERTX = V14*V23 -V13*V24
-C
-         VERTEX = DCMPLX( DVERTX ) * (G*G)
-C
-      RETURN
-      END
-C
-C ======================================================================
-C
-      SUBROUTINE GGGXXX(WM,WP,W3,G , VERTEX)
-C
-C this subroutine computes an amplitude of the three-point coupling of  
-C the gauge bosons.                                                     
-C                                                                       
-C input:                                                                
-C       complex wm(6)          : vector               flow-out w-       
-C       complex wp(6)          : vector               flow-out w+       
-C       complex w3(6)          : vector               j3 or a    or z   
-C       real    g              : coupling constant    gw or gwwa or gwwz
-C                                                                       
-C output:                                                               
-C       complex vertex         : amplitude               gamma(wm,wp,w3)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 WM(6),WP(6),W3(6),VERTEX, 
-     &        XV1,XV2,XV3,V12,V23,V31,P12,P13,P21,P23,P31,P32
-      REAL*8    PWM(0:3),PWP(0:3),PW3(0:3),G
-      REAL*8 RXZERO, RTENTH
-      PARAMETER( RXZERO=0.0D0, RTENTH=0.1D0 )
-C
-      PWM(0)=DBLE( WM(5))
-      PWM(1)=DBLE( WM(6))
-      PWM(2)=DIMAG(WM(6))
-      PWM(3)=DIMAG(WM(5))
-      PWP(0)=DBLE( WP(5))
-      PWP(1)=DBLE( WP(6))
-      PWP(2)=DIMAG(WP(6))
-      PWP(3)=DIMAG(WP(5))
-      PW3(0)=DBLE( W3(5))
-      PW3(1)=DBLE( W3(6))
-      PW3(2)=DIMAG(W3(6))
-      PW3(3)=DIMAG(W3(5))
-C
-      V12=WM(1)*WP(1)-WM(2)*WP(2)-WM(3)*WP(3)-WM(4)*WP(4)
-      V23=WP(1)*W3(1)-WP(2)*W3(2)-WP(3)*W3(3)-WP(4)*W3(4)
-      V31=W3(1)*WM(1)-W3(2)*WM(2)-W3(3)*WM(3)-W3(4)*WM(4)
-      XV1=RXZERO
-      XV2=RXZERO
-      XV3=RXZERO
-      IF ( ABS(WM(1)) .NE. RXZERO ) THEN
-         IF (ABS(WM(1)).GE.MAX(ABS(WM(2)),ABS(WM(3)),ABS(WM(4)))
-     $        *RTENTH)
-     &      XV1=PWM(0)/WM(1)
-      ENDIF
-      IF ( ABS(WP(1)) .NE. RXZERO) THEN
-         IF (ABS(WP(1)).GE.MAX(ABS(WP(2)),ABS(WP(3)),ABS(WP(4)))
-     $        *RTENTH)
-     &      XV2=PWP(0)/WP(1)
-      ENDIF
-      IF ( ABS(W3(1)) .NE. RXZERO) THEN
-         IF ( ABS(W3(1)).GE.MAX(ABS(W3(2)),ABS(W3(3)),ABS(W3(4)))
-     $        *RTENTH)
-     &      XV3=PW3(0)/W3(1)
-      ENDIF
-      P12= (PWM(0)-XV1*WM(1))*WP(1)-(PWM(1)-XV1*WM(2))*WP(2)
-     &    -(PWM(2)-XV1*WM(3))*WP(3)-(PWM(3)-XV1*WM(4))*WP(4)
-      P13= (PWM(0)-XV1*WM(1))*W3(1)-(PWM(1)-XV1*WM(2))*W3(2)
-     &    -(PWM(2)-XV1*WM(3))*W3(3)-(PWM(3)-XV1*WM(4))*W3(4)
-      P21= (PWP(0)-XV2*WP(1))*WM(1)-(PWP(1)-XV2*WP(2))*WM(2)
-     &    -(PWP(2)-XV2*WP(3))*WM(3)-(PWP(3)-XV2*WP(4))*WM(4)
-      P23= (PWP(0)-XV2*WP(1))*W3(1)-(PWP(1)-XV2*WP(2))*W3(2)
-     &    -(PWP(2)-XV2*WP(3))*W3(3)-(PWP(3)-XV2*WP(4))*W3(4)
-      P31= (PW3(0)-XV3*W3(1))*WM(1)-(PW3(1)-XV3*W3(2))*WM(2)
-     &    -(PW3(2)-XV3*W3(3))*WM(3)-(PW3(3)-XV3*W3(4))*WM(4)
-      P32= (PW3(0)-XV3*W3(1))*WP(1)-(PW3(1)-XV3*W3(2))*WP(2)
-     &    -(PW3(2)-XV3*W3(3))*WP(3)-(PW3(3)-XV3*W3(4))*WP(4)
-C
-      VERTEX = -(V12*(P13-P23)+V23*(P21-P31)+V31*(P32-P12))*G
-C
-      RETURN
-      END
-      SUBROUTINE HIOXXX(FI,FO,GC,SMASS,SWIDTH , HIO)
-C
-C this subroutine computes an off-shell scalar current from an external
-C fermion pair.
-C       
-C input:
-C       complex fi(6)          : flow-in  fermion                   |fi>
-C       complex fo(6)          : flow-out fermion                   <fo|
-C       complex gc(2)          : coupling constants                 gchf
-C       real    smass          : mass  of output scalar s
-C       real    swidth         : width of output scalar s
-C       
-C output:
-C       complex hio(3)         : scalar current             j(<fi|s|fo>)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FI(6),FO(6),HIO(3),GC(2),DN
-      REAL*8  Q(0:3),SMASS,SWIDTH,Q2
-C       
-      HIO(2) = FO(5)-FI(5)
-      HIO(3) = FO(6)-FI(6)
-C       
-      Q(0)=DBLE( HIO(2))
-      Q(1)=DBLE( HIO(3))
-      Q(2)=DIMAG(HIO(3))
-      Q(3)=DIMAG(HIO(2))
-      Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-C
-      DN=-DCMPLX(Q2-SMASS**2,DMAX1(DSIGN(SMASS*SWIDTH,Q2),0.D0))
-C
-      HIO(1) = ( GC(1)*(FO(1)*FI(1)+FO(2)*FI(2))
-     &          +GC(2)*(FO(3)*FI(3)+FO(4)*FI(4)) )/DN
-C
-      RETURN
-      END
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE HSSSXX(S1,S2,S3,G,SMASS,SWIDTH , HSSS)
-C
-C This subroutine computes an off-shell scalar current from the four-   
-C scalar coupling.                                                      
-C                                                                       
-C INPUT:                                                                
-C       complex S1(3)          : first  scalar                        S1
-C       complex S2(3)          : second scalar                        S2
-C       complex S3(3)          : third  scalar                        S3
-C       real    G              : coupling constant                 GHHHH
-C       real    SMASS          : mass  of OUTPUT scalar S'              
-C       real    SWIDTH         : width of OUTPUT scalar S'              
-C                                                                       
-C OUTPUT:                                                               
-C       complex HSSS(3)        : scalar current           J(S':S1,S2,S3)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 S1(3),S2(3),S3(3),HSSS(3),DG
-      REAL*8     Q(0:3),G,SMASS,SWIDTH,Q2
-C
-      HSSS(2) = S1(2)+S2(2)+S3(2)
-      HSSS(3) = S1(3)+S2(3)+S3(3)
-C
-      Q(0)=DBLE( HSSS(2))
-      Q(1)=DBLE( HSSS(3))
-      Q(2)=DIMAG(HSSS(3))
-      Q(3)=DIMAG(HSSS(2))
-      Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-C
-      DG=-G/DCMPLX( Q2-SMASS**2,MAX(SIGN(SMASS*SWIDTH ,Q2),0.D0))
-C
-      HSSS(1) = DG * S1(1)*S2(1)*S3(1)
-C
-      RETURN
-      END
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE HSSXXX(S1,S2,G,SMASS,SWIDTH , HSS)
-C
-C This subroutine computes an off-shell scalar current from the three-  
-C scalar coupling.                                                      
-C                                                                       
-C INPUT:                                                                
-C       complex S1(3)          : first  scalar                        S1
-C       complex S2(3)          : second scalar                        S2
-C       real    G              : coupling constant                  GHHH
-C       real    SMASS          : mass  of OUTPUT scalar S'              
-C       real    SWIDTH         : width of OUTPUT scalar S'              
-C                                                                       
-C OUTPUT:                                                               
-C       complex HSS(3)         : scalar current              J(S':S1,S2)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 S1(3),S2(3),HSS(3),DG
-      REAL*8  Q(0:3),G,SMASS,SWIDTH,Q2
-C
-      HSS(2) = S1(2)+S2(2)
-      HSS(3) = S1(3)+S2(3)
-C
-      Q(0)=DBLE( HSS(2))
-      Q(1)=DBLE( HSS(3))
-      Q(2)=DIMAG(HSS(3))
-      Q(3)=DIMAG(HSS(2))
-      Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-C
-      DG=-G/DCMPLX( Q2-SMASS**2, MAX(SIGN(SMASS*SWIDTH ,Q2),0.D0))
-C
-      HSS(1) = DG*S1(1)*S2(1)
-C
-      RETURN
-      END
-C
-C ======================================================================
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE HVSXXX(VC,SC,G,SMASS,SWIDTH , HVS)
-C
-C this subroutine computes an off-shell scalar current from the vector- 
-C scalar-scalar coupling.  the coupling is absent in the minimal sm in  
-C unitary gauge.                                                        
-C                                                                       
-C input:                                                                
-C       complex vc(6)          : input vector                          v
-C       complex sc(3)          : input scalar                          s
-C       complex g              : coupling constant (s charge)           
-C       real    smass          : mass  of output scalar s'              
-C       real    swidth         : width of output scalar s'              
-C                                                                       
-C examples of the coupling constant g for susy particles are as follows:
-C   -----------------------------------------------------------         
-C   |    s1    | (q,i3) of s1  ||   v=a   |   v=z   |   v=w   |         
-C   -----------------------------------------------------------         
-C   | nu~_l    | (  0  , +1/2) ||   ---   |  gzn(1) |  gwf(1) |         
-C   | e~_l     | ( -1  , -1/2) ||  gal(1) |  gzl(1) |  gwf(1) |         
-C   | u~_l     | (+2/3 , +1/2) ||  gau(1) |  gzu(1) |  gwf(1) |         
-C   | d~_l     | (-1/3 , -1/2) ||  gad(1) |  gzd(1) |  gwf(1) |         
-C   -----------------------------------------------------------         
-C   | e~_r-bar | ( +1  ,  0  ) || -gal(2) | -gzl(2) | -gwf(2) |         
-C   | u~_r-bar | (-2/3 ,  0  ) || -gau(2) | -gzu(2) | -gwf(2) |         
-C   | d~_r-bar | (+1/3 ,  0  ) || -gad(2) | -gzd(2) | -gwf(2) |         
-C   -----------------------------------------------------------         
-C where the sc charge is defined by the flowing-out quantum number.     
-C                                                                       
-C output:                                                               
-C       complex hvs(3)         : scalar current                j(s':v,s)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 VC(6),SC(3),HVS(3),DG,QVV,QPV,G
-      REAL*8    QV(0:3),QP(0:3),QA(0:3),SMASS,SWIDTH,Q2
-C
-      HVS(2) = VC(5)+SC(2)
-      HVS(3) = VC(6)+SC(3)
-C
-      QV(0)=DBLE(  VC(5))
-      QV(1)=DBLE(  VC(6))
-      QV(2)=DIMAG( VC(6))
-      QV(3)=DIMAG( VC(5))
-      QP(0)=DBLE(  SC(2))
-      QP(1)=DBLE(  SC(3))
-      QP(2)=DIMAG( SC(3))
-      QP(3)=DIMAG( SC(2))
-      QA(0)=DBLE( HVS(2))
-      QA(1)=DBLE( HVS(3))
-      QA(2)=DIMAG(HVS(3))
-      QA(3)=DIMAG(HVS(2))
-      Q2=QA(0)**2-(QA(1)**2+QA(2)**2+QA(3)**2)
-C
-      DG=-G/DCMPLX( Q2-SMASS**2 , MAX(DSIGN( SMASS*SWIDTH ,Q2),0D0) )
-      QVV=QV(0)*VC(1)-QV(1)*VC(2)-QV(2)*VC(3)-QV(3)*VC(4)
-      QPV=QP(0)*VC(1)-QP(1)*VC(2)-QP(2)*VC(3)-QP(3)*VC(4)
-C
-      HVS(1) = DG*(2D0*QPV+QVV)*SC(1)
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE HVVXXX(V1,V2,G,SMASS,SWIDTH , HVV)
-C
-C this subroutine computes an off-shell scalar current from the vector- 
-C vector-scalar coupling.                                               
-C                                                                       
-C input:                                                                
-C       complex v1(6)          : first  vector                        v1
-C       complex v2(6)          : second vector                        v2
-C       real    g              : coupling constant                  gvvh
-C       real    smass          : mass  of output scalar s               
-C       real    swidth         : width of output scalar s               
-C                                                                       
-C output:                                                               
-C       complex hvv(3)         : off-shell scalar current     j(s:v1,v2)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 V1(6),V2(6),HVV(3),DG
-      REAL*8    Q(0:3),G,SMASS,SWIDTH,Q2
-      REAL*8 RXZERO
-      PARAMETER( RXZERO=0.0D0 )
-C
-      HVV(2) = V1(5)+V2(5)
-      HVV(3) = V1(6)+V2(6)
-C
-      Q(0)=DBLE( HVV(2))
-      Q(1)=DBLE( HVV(3))
-      Q(2)=DIMAG(HVV(3))
-      Q(3)=DIMAG(HVV(2))
-      Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-C
-      DG=-G/DCMPLX( Q2-SMASS**2 , MAX(SIGN( SMASS*SWIDTH ,Q2),RXZERO) )
-C
-      HVV(1) = DG*(V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4))
-C
-      RETURN
-      END
-C
-C ======================================================================
-C
-      SUBROUTINE IOSXXX(FI,FO,SC,GC , VERTEX)
-C
-C This subroutine computes an amplitude of the fermion-fermion-scalar   
-C coupling.                                                             
-C                                                                       
-C INPUT:                                                                
-C       complex FI(6)          : flow-in  fermion                   |FI>
-C       complex FO(6)          : flow-out fermion                   <FO|
-C       complex SC(3)          : input    scalar                      S 
-C       complex GC(2)          : coupling constants                 GCHF
-C                                                                       
-C OUTPUT:                                                               
-C       complex VERTEX         : amplitude                     <FO|S|FI>
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FI(6),FO(6),SC(3),GC(2),VERTEX
-C
-      VERTEX = SC(1)*( GC(1)*(FI(1)*FO(1)+FI(2)*FO(2))
-     &                +GC(2)*(FI(3)*FO(3)+FI(4)*FO(4)) )
-C
-      RETURN          
-      END
-C
-C ======================================================================
-C
-      SUBROUTINE IOVXXX(FI,FO,VC,G , VERTEX)
-C
-C this subroutine computes an amplitude of the fermion-fermion-vector   
-C coupling.                                                             
-C                                                                       
-C input:                                                                
-C       complex fi(6)          : flow-in  fermion                   |fi>
-C       complex fo(6)          : flow-out fermion                   <fo|
-C       complex vc(6)          : input    vector                      v 
-C       real    g(2)           : coupling constants                  gvf
-C                                                                       
-C output:                                                               
-C       complex vertex         : amplitude                     <fo|v|fi>
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FI(6),FO(6),VC(6),VERTEX
-      REAL*8    G(2)
-      REAL*8 RXZERO, RXONE
-      PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
-      COMPLEX*16 CXIMAG
-      LOGICAL FIRST
-      SAVE CXIMAG,FIRST
-      DATA FIRST/.TRUE./
-C
-C          Fix compilation with g77
-      IF(FIRST) THEN
-        FIRST=.FALSE.
-        CXIMAG=DCMPLX( RXZERO, RXONE )
-      ENDIF
-C
-
-      VERTEX =  G(1)*( (FO(3)*FI(1)+FO(4)*FI(2))*VC(1)
-     &                +(FO(3)*FI(2)+FO(4)*FI(1))*VC(2)
-     &                -(FO(3)*FI(2)-FO(4)*FI(1))*VC(3)*CXIMAG
-     &                +(FO(3)*FI(1)-FO(4)*FI(2))*VC(4)        )
-C
-      IF ( G(2) .NE. RXZERO ) THEN
-         VERTEX = VERTEX
-     &          + G(2)*( (FO(1)*FI(3)+FO(2)*FI(4))*VC(1)
-     &                  -(FO(1)*FI(4)+FO(2)*FI(3))*VC(2)
-     &                  +(FO(1)*FI(4)-FO(2)*FI(3))*VC(3)*CXIMAG
-     &                  -(FO(1)*FI(3)-FO(2)*FI(4))*VC(4)        )
-      END IF
-C
-      RETURN
-      END
-C
-C      Subroutine returns the desired fermion or
-C      anti-fermion spinor. ie., |f>
-C      A replacement for the HELAS routine IXXXXX
-C
-C      Adam Duff,  1992 August 31
-C      <duff@phenom.physics.wisc.edu>
-C
-      SUBROUTINE IXXXXX(P,FMASS,NHEL,NSF,FI)
-C          P        IN: FOUR VECTOR MOMENTUM
-C          FMASS    IN: FERMION MASS
-C          NHEL     IN: SPINOR HELICITY, -1 OR 1
-C          NSF      IN: -1=ANTIFERMION, 1=FERMION
-C          FI       OUT: FERMION WAVEFUNCTION
-C
-C declare input/output variables
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FI(6)
-      INTEGER*4 NHEL, NSF
-      REAL*8 P(0:3), FMASS
-      REAL*8 RXZERO, RXONE, RXTWO
-      PARAMETER( RXZERO=0.0D0, RXONE=1.0D0, RXTWO=2.0D0 )
-      REAL*8 PLAT, PABS, OMEGAP, OMEGAM, RS2PA, SPAZ
-      COMPLEX*16 CXZERO
-C
-C declare local variables
-C
-      LOGICAL FIRST
-      SAVE CXZERO,FIRST
-      DATA FIRST/.TRUE./
-C
-C          Fix compilation with g77
-      IF(FIRST) THEN
-        FIRST=.FALSE.
-        CXZERO=DCMPLX( RXZERO, RXZERO )
-      ENDIF
-C
-C define kinematic parameters
-C
-      FI(5) = DCMPLX( P(0), P(3) ) * NSF
-      FI(6) = DCMPLX( P(1), P(2) ) * NSF
-      PLAT = SQRT( P(1)**2 + P(2)**2 )
-      PABS = SQRT( P(1)**2 + P(2)**2 + P(3)**2 )
-      OMEGAP = SQRT( P(0) + PABS )
-C
-C do massive fermion case
-C
-      IF ( FMASS .NE. RXZERO ) THEN
-         OMEGAM = FMASS / OMEGAP
-         IF ( NSF .EQ. 1 ) THEN
-            IF ( NHEL .EQ. 1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = DCMPLX( OMEGAM, RXZERO )
-                     FI(2) = CXZERO
-                     FI(3) = DCMPLX( OMEGAP, RXZERO )
-                     FI(4) = CXZERO
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS + P(3) )
-                     FI(1) = OMEGAM * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FI(2) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( P(1), P(2) )
-                     FI(3) = OMEGAP * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FI(4) = OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( P(1), P(2) )
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = CXZERO
-                     FI(2) = DCMPLX( OMEGAM, RXZERO )
-                     FI(3) = CXZERO
-                     FI(4) = DCMPLX( OMEGAP, RXZERO )
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS - P(3) )
-                     FI(1) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FI(2) = OMEGAM * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( P(1), P(2) )
-                     FI(3) = OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FI(4) = OMEGAP * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( P(1), P(2) )
-                  END IF
-               END IF
-            ELSE IF ( NHEL .EQ. -1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = CXZERO
-                     FI(2) = DCMPLX( OMEGAP, RXZERO )
-                     FI(3) = CXZERO
-                     FI(4) = DCMPLX( OMEGAM, RXZERO )
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS + P(3) )
-                     FI(1) = OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(2) = OMEGAP * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FI(3) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(4) = OMEGAM * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = DCMPLX( -OMEGAP, RXZERO )
-                     FI(2) = CXZERO
-                     FI(3) = DCMPLX( -OMEGAM, RXZERO )
-                     FI(4) = CXZERO
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS - P(3) )
-                     FI(1) = OMEGAP * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(2) = OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FI(3) = OMEGAM * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(4) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                  END IF
-               END IF
-            ELSE
-               STOP 'IXXXXX:  FERMION HELICITY MUST BE +1,-1'
-            END IF
-         ELSE IF ( NSF .EQ. -1 ) THEN
-            IF ( NHEL .EQ. 1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = CXZERO
-                     FI(2) = DCMPLX( -OMEGAP, RXZERO )
-                     FI(3) = CXZERO
-                     FI(4) = DCMPLX( OMEGAM, RXZERO )
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS + P(3) )
-                     FI(1) = -OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(2) = -OMEGAP * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FI(3) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(4) = OMEGAM * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = DCMPLX( OMEGAP, RXZERO )
-                     FI(2) = CXZERO
-                     FI(3) = DCMPLX( -OMEGAM, RXZERO )
-                     FI(4) = CXZERO
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS - P(3) )
-                     FI(1) = -OMEGAP * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(2) = -OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FI(3) = OMEGAM * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(4) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                  END IF
-               END IF
-            ELSE IF ( NHEL .EQ. -1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = DCMPLX( OMEGAM, RXZERO )
-                     FI(2) = CXZERO
-                     FI(3) = DCMPLX( -OMEGAP, RXZERO )
-                     FI(4) = CXZERO
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS + P(3) )
-                     FI(1) = OMEGAM * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FI(2) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( P(1), P(2) )
-                     FI(3) = -OMEGAP * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FI(4) = -OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( P(1), P(2) )
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = CXZERO
-                     FI(2) = DCMPLX( OMEGAM, RXZERO )
-                     FI(3) = CXZERO
-                     FI(4) = DCMPLX( -OMEGAP, RXZERO )
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS - P(3) )
-                     FI(1) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FI(2) = OMEGAM * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( P(1), P(2) )
-                     FI(3) = -OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FI(4) = -OMEGAP * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( P(1), P(2) )
-                  END IF
-               END IF
-            ELSE
-               STOP 'IXXXXX:  FERMION HELICITY MUST BE +1,-1'
-            END IF
-         ELSE
-            STOP 'IXXXXX:  FERMION TYPE MUST BE +1,-1'
-         END IF
-C
-C do massless fermion case
-C
-      ELSE
-         IF ( NSF .EQ. 1 ) THEN
-            IF ( NHEL .EQ. 1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = CXZERO
-                     FI(2) = CXZERO
-                     FI(3) = DCMPLX( OMEGAP, RXZERO )
-                     FI(4) = CXZERO
-                  ELSE
-                     SPAZ = SQRT( PABS + P(3) )
-                     FI(1) = CXZERO
-                     FI(2) = CXZERO
-                     FI(3) = DCMPLX( SPAZ, RXZERO )
-                     FI(4) = RXONE / SPAZ
-     &                     * DCMPLX( P(1), P(2) )
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = CXZERO
-                     FI(2) = CXZERO
-                     FI(3) = CXZERO
-                     FI(4) = DCMPLX( OMEGAP, RXZERO )
-                  ELSE
-                     SPAZ = SQRT( PABS - P(3) )
-                     FI(1) = CXZERO
-                     FI(2) = CXZERO
-                     FI(3) = RXONE / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FI(4) = SPAZ / PLAT
-     &                     * DCMPLX( P(1), P(2) )
-                  END IF
-               END IF
-            ELSE IF ( NHEL .EQ. -1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = CXZERO
-                     FI(2) = DCMPLX( OMEGAP, RXZERO )
-                     FI(3) = CXZERO
-                     FI(4) = CXZERO
-                  ELSE
-                     SPAZ = SQRT( PABS + P(3) )
-                     FI(1) = RXONE / SPAZ
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(2) = DCMPLX( SPAZ, RXZERO )
-                     FI(3) = CXZERO
-                     FI(4) = CXZERO
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = DCMPLX( -OMEGAP, RXZERO )
-                     FI(2) = CXZERO
-                     FI(3) = CXZERO
-                     FI(4) = CXZERO
-                  ELSE
-                     SPAZ = SQRT( PABS - P(3) )
-                     FI(1) = SPAZ / PLAT
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(2) = RXONE / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FI(3) = CXZERO
-                     FI(4) = CXZERO
-                  END IF
-               END IF
-            ELSE
-               STOP 'IXXXXX:  FERMION HELICITY MUST BE +1,-1'
-            END IF
-         ELSE IF ( NSF .EQ. -1 ) THEN
-            IF ( NHEL .EQ. 1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = CXZERO
-                     FI(2) = DCMPLX( -OMEGAP, RXZERO )
-                     FI(3) = CXZERO
-                     FI(4) = CXZERO
-                  ELSE
-                     SPAZ = SQRT( PABS + P(3) )
-                     FI(1) = -RXONE / SPAZ
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(2) = DCMPLX( -SPAZ, RXZERO )
-                     FI(3) = CXZERO
-                     FI(4) = CXZERO
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = DCMPLX( OMEGAP, RXZERO )
-                     FI(2) = CXZERO
-                     FI(3) = CXZERO
-                     FI(4) = CXZERO
-                  ELSE
-                     SPAZ = SQRT( PABS - P(3) )
-                     FI(1) = -SPAZ / PLAT
-     &                     * DCMPLX( -P(1), P(2) )
-                     FI(2) = -RXONE / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FI(3) = CXZERO
-                     FI(4) = CXZERO
-                  END IF
-               END IF
-            ELSE IF ( NHEL .EQ. -1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = CXZERO
-                     FI(2) = CXZERO
-                     FI(3) = DCMPLX( -OMEGAP, RXZERO )
-                     FI(4) = CXZERO
-                  ELSE
-                     SPAZ = SQRT( PABS + P(3) )
-                     FI(1) = CXZERO
-                     FI(2) = CXZERO
-                     FI(3) = DCMPLX( -SPAZ, RXZERO )
-                     FI(4) = -RXONE / SPAZ
-     &                     * DCMPLX( P(1), P(2) )
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FI(1) = CXZERO
-                     FI(2) = CXZERO
-                     FI(3) = CXZERO
-                     FI(4) = DCMPLX( -OMEGAP, RXZERO )
-                  ELSE
-                     SPAZ = SQRT( PABS - P(3) )
-                     FI(1) = CXZERO
-                     FI(2) = CXZERO
-                     FI(3) = -RXONE / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FI(4) = -SPAZ / PLAT
-     &                     * DCMPLX( P(1), P(2) )
-                  END IF
-               END IF
-            ELSE
-               STOP 'IXXXXX:  FERMION HELICITY MUST BE +1,-1'
-            END IF
-         ELSE
-            STOP 'IXXXXX:  FERMION TYPE MUST BE +1,-1'
-         END IF
-      END IF
-C
-C done
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE J3XXXX(FI,FO,GAF,GZF,ZMASS,ZWIDTH , J3)
-C
-C this subroutine computes the sum of photon and z currents with the    
-C suitable weights ( j(w3) = cos(theta_w) j(z) + sin(theta_w) j(a) ).   
-C the output j3 is useful as an input of vvvxxx, jvvxxx or w3w3xx.      
-C the photon propagator is given in feynman gauge, and the z propagator 
-C is given in unitary gauge.                                            
-C                                                                       
-C input:                                                                
-C       complex fi(6)          : flow-in  fermion                   |fi>
-C       complex fo(6)          : flow-out fermion                   <fo|
-C       real    gaf(2)         : fi couplings with a                 gaf
-C       real    gzf(2)         : fi couplings with z                 gzf
-C       real    zmass          : mass  of z                             
-C       real    zwidth         : width of z                             
-C                                                                       
-C output:                                                               
-C       complex j3(6)          : w3 current             j^mu(<fo|w3|fi>)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FI(6),FO(6),J3(6),
-     &        C0L,C1L,C2L,C3L,CSL,C0R,C1R,C2R,C3R,CSR,DZ,DDIF
-      REAL*8    GAF(2),GZF(2),Q(0:3),ZMASS,ZWIDTH,ZM2,ZMW,Q2,DA,WW,
-     &        CW,SW,GN,GZ3L,GA3L
-C
-      REAL*8 RXZERO, RXONE
-      PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
-      COMPLEX*16 CXIMAG
-      LOGICAL FIRST
-      SAVE CXIMAG,FIRST
-      DATA FIRST/.TRUE./
-C
-C          Fix compilation with g77
-      IF(FIRST) THEN
-        FIRST=.FALSE.
-        CXIMAG=DCMPLX( RXZERO, RXONE )
-      ENDIF
-C
-      J3(5) = FO(5)-FI(5)
-      J3(6) = FO(6)-FI(6)
-C
-      Q(0)=-DBLE( J3(5))
-      Q(1)=-DBLE( J3(6))
-      Q(2)=-DIMAG(J3(6))
-      Q(3)=-DIMAG(J3(5))
-      Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-      ZM2=ZMASS**2
-      ZMW=ZMASS*ZWIDTH
-C
-      DA=RXONE/Q2
-      WW=MAX(DSIGN( ZMW ,Q2),RXZERO)
-      DZ=RXONE/DCMPLX( Q2-ZM2 , WW )
-      DDIF=DCMPLX( -ZM2 , WW )*DA*DZ
-C
-C ddif is the difference : ddif=da-dz
-C  for the running width, use below instead of the above ww,dz and ddif.
-C      ww=max( zwidth*q2/zmass ,r_zero)
-C      dz=r_one/dcmplx( q2-zm2 , ww )
-C      ddif=dcmplx( -zm2 , ww )*da*dz
-C
-      CW=RXONE/SQRT(RXONE+(GZF(2)/GAF(2))**2)
-      SW=SQRT((RXONE-CW)*(RXONE+CW))
-      GN=GAF(2)*SW
-      GZ3L=GZF(1)*CW
-      GA3L=GAF(1)*SW
-      C0L=  FO(3)*FI(1)+FO(4)*FI(2)
-      C0R=  FO(1)*FI(3)+FO(2)*FI(4)
-      C1L=-(FO(3)*FI(2)+FO(4)*FI(1))
-      C1R=  FO(1)*FI(4)+FO(2)*FI(3)
-      C2L= (FO(3)*FI(2)-FO(4)*FI(1))*CXIMAG
-      C2R=(-FO(1)*FI(4)+FO(2)*FI(3))*CXIMAG
-      C3L= -FO(3)*FI(1)+FO(4)*FI(2)
-      C3R=  FO(1)*FI(3)-FO(2)*FI(4)
-      CSL=(Q(0)*C0L-Q(1)*C1L-Q(2)*C2L-Q(3)*C3L)/ZM2
-      CSR=(Q(0)*C0R-Q(1)*C1R-Q(2)*C2R-Q(3)*C3R)/ZM2
-C
-      J3(1) =  GZ3L*DZ*(C0L-CSL*Q(0))+GA3L*C0L*DA
-     &       + GN*(C0R*DDIF-CSR*Q(0)*DZ)
-      J3(2) =  GZ3L*DZ*(C1L-CSL*Q(1))+GA3L*C1L*DA
-     &       + GN*(C1R*DDIF-CSR*Q(1)*DZ)
-      J3(3) =  GZ3L*DZ*(C2L-CSL*Q(2))+GA3L*C2L*DA
-     &       + GN*(C2R*DDIF-CSR*Q(2)*DZ)
-      J3(4) =  GZ3L*DZ*(C3L-CSL*Q(3))+GA3L*C3L*DA
-     &       + GN*(C3R*DDIF-CSR*Q(3)*DZ)
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JEEXXX(EB,EF,SHLF,CHLF,PHI,NHB,NHF,NSF , JEE)
-C
-C This subroutine computes an off-shell photon wavefunction emitted from
-C the electron or positron beam, with a special care for the small angle
-C region.  The momenta are measured in the laboratory frame, where the  
-C e- (e+) beam is along the positive (negative) z axis.                 
-C                                                                       
-C INPUT:                                                                
-C       real    EB             : energy (GeV)    of beam  e-/e+         
-C       real    EF             : energy (GeV)    of final e-/e+         
-C       real    SHLF           : sin(theta/2)    of final e-/e+         
-C       real    CHLF           : cos(theta/2)    of final e-/e+         
-C       real    PHI            : azimuthal angle of final e-/e+         
-C       integer NHB  = -1 or 1 : helicity        of beam  e-/e+         
-C       integer NHF  = -1 or 1 : helicity        of final e-/e+         
-C       integer NSF  = -1 or 1 : +1 for electron, -1 for positron       
-C                                                                       
-C OUTPUT:                                                               
-C       complex JEE(6)         : off-shell photon          J^mu(<e|A|e>)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 JEE(6),COEFF
-      REAL*8  CS(2),EB,EF,SHLF,CHLF,PHI,ME,ALPHA,GAL,HI,SF,SFH,X,ME2,Q2,
-     &        RFP,RFM,SNP,CSP,RXC,C,S
-      INTEGER NHB,NHF,NSF
-C
-      ME   =0.51099906D-3
-      ALPHA=1./128.
-      GAL  =SQRT(ALPHA*4.*3.14159265D0)
-C
-      HI =NHB
-      SF =NSF
-      SFH=NHB*NSF
-      CS((3+NSF)/2)=SHLF
-      CS((3-NSF)/2)=CHLF
-C CS(1)=CHLF and CS(2)=SHLF for electron
-C CS(1)=SHLF and CS(2)=CHLF for positron
-      X=EF/EB
-      ME2=ME**2
-      Q2=-4.*CS(2)**2*(EF*EB-ME2)
-     &   +SF*(1.-X)**2/X*(SHLF+CHLF)*(SHLF-CHLF)*ME2
-      RFP=(1+NSF)
-      RFM=(1-NSF)
-      SNP=SIN(PHI)
-      CSP=COS(PHI)
-C
-      IF (NHB.EQ.NHF) THEN
-         RXC=2.*X/(1.-X)*CS(1)**2
-         COEFF= GAL*2.*EB*SQRT(X)*CS(2)/Q2
-     &         *(DCMPLX( RFP )-RFM*DCMPLX( CSP ,-SNP*HI ))*.5
-         JEE(1) =  DCMPLX( 0.D0 )
-         JEE(2) =  COEFF*DCMPLX( (1.+RXC)*CSP ,-SFH*SNP )
-         JEE(3) =  COEFF*DCMPLX( (1.+RXC)*SNP , SFH*CSP )
-         JEE(4) =  COEFF*(-SF*RXC/CS(1)*CS(2))
-      ELSE
-         COEFF= GAL*ME/Q2/SQRT(X)
-     &         *(DCMPLX( RFP )+RFM*DCMPLX( CSP , SNP*HI ))*.5*HI
-         JEE(1) = -COEFF*(1.+X)*CS(2)*DCMPLX( CSP , SFH*SNP )
-         JEE(2) =  COEFF*(1.-X)*CS(1)
-         JEE(3) =  JEE(2)*DCMPLX( 0.D0 , SFH )
-         JEE(4) =  JEE(1)*SF*(1.-X)/(1.+X)
-      ENDIF
-C
-      C=(CHLF+SHLF)*(CHLF-SHLF)
-      S=2.*CHLF*SHLF
-C
-      JEE(5) = -EB*DCMPLX( 1.-X , SF-X*C )
-      JEE(6) =  EB*X*S*DCMPLX( CSP , SNP )
-C
-      RETURN          
-      END
-C
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JGGGXX(W1,W2,W3,G, JW3W)
-C
-C this subroutine computes an off-shell w+, w-, w3, z or photon current 
-C from the four-point gauge boson coupling, including the contributions 
-C of w exchange diagrams.  the vector propagator is given in feynman    
-C gauge for a photon and in unitary gauge for w and z bosons.  if one   
-C sets wmass=0.0, then the ggg-->g current is given (see sect 2.9.1 of 
-C the manual).                                                          
-C                                                                       
-C input:                                                                
-C       complex w1(6)          : first  vector                        w1
-C       complex w2(6)          : second vector                        w2
-C       complex w3(6)          : third  vector                        w3
-C       real    g             : first  coupling constant               
-C                                                  (see the table below)
-C                                                                       
-C output:                                                               
-C       complex jw3w(6)        : w current             j^mu(w':w1,w2,w3)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16  W1(6),W2(6),W3(6),JW3W(6)
-      COMPLEX*16 DW1(0:3),DW2(0:3),DW3(0:3),
-     &           JJ(0:3),DV,W32,W13
-      REAL*8     P1(0:3),P2(0:3),P3(0:3),Q(0:3),G,DG2,Q2
-C
-      REAL*8 RXZERO
-      PARAMETER( RXZERO=0.0D0 )
-C
-      JW3W(5) = W1(5)+W2(5)+W3(5)
-      JW3W(6) = W1(6)+W2(6)+W3(6)
-C
-      DW1(0)=DCMPLX(W1(1))
-      DW1(1)=DCMPLX(W1(2))
-      DW1(2)=DCMPLX(W1(3))
-      DW1(3)=DCMPLX(W1(4))
-      DW2(0)=DCMPLX(W2(1))
-      DW2(1)=DCMPLX(W2(2))
-      DW2(2)=DCMPLX(W2(3))
-      DW2(3)=DCMPLX(W2(4))
-      DW3(0)=DCMPLX(W3(1))
-      DW3(1)=DCMPLX(W3(2))
-      DW3(2)=DCMPLX(W3(3))
-      DW3(3)=DCMPLX(W3(4))
-      P1(0)=DBLE(      W1(5))
-      P1(1)=DBLE(      W1(6))
-      P1(2)=DBLE(DIMAG(W1(6)))
-      P1(3)=DBLE(DIMAG(W1(5)))
-      P2(0)=DBLE(      W2(5))
-      P2(1)=DBLE(      W2(6))
-      P2(2)=DBLE(DIMAG(W2(6)))
-      P2(3)=DBLE(DIMAG(W2(5)))
-      P3(0)=DBLE(      W3(5))
-      P3(1)=DBLE(      W3(6))
-      P3(2)=DBLE(DIMAG(W3(6)))
-      P3(3)=DBLE(DIMAG(W3(5)))
-      Q(0)=-(P1(0)+P2(0)+P3(0))
-      Q(1)=-(P1(1)+P2(1)+P3(1))
-      Q(2)=-(P1(2)+P2(2)+P3(2))
-      Q(3)=-(P1(3)+P2(3)+P3(3))
-
-      Q2 =Q(0)**2 -(Q(1)**2 +Q(2)**2 +Q(3)**2)
-
-      DG2=DBLE(G)*DBLE(G)
-C
-      DV = 1.0D0/DCMPLX( Q2 )
-
-C  for the running width, use below instead of the above dv.
-C      dv = 1.0d0/dcmplx( q2 -mv2 , dmax1(dwv*q2/dmv,0.d0) )
-C
-      W32=DW3(0)*DW2(0)-DW3(1)*DW2(1)-DW3(2)*DW2(2)-DW3(3)*DW2(3)
-C
-C     
-      W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
-C     
-      JJ(0)=DG2*( DW1(0)*W32 - DW2(0)*W13 )
-      JJ(1)=DG2*( DW1(1)*W32 - DW2(1)*W13 )
-      JJ(2)=DG2*( DW1(2)*W32 - DW2(2)*W13 )
-      JJ(3)=DG2*( DW1(3)*W32 - DW2(3)*W13 )
-C     
-      JW3W(1) = DCMPLX( JJ(0)*DV )
-      JW3W(2) = DCMPLX( JJ(1)*DV )
-      JW3W(3) = DCMPLX( JJ(2)*DV )
-      JW3W(4) = DCMPLX( JJ(3)*DV )
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JGGXXX(V1,V2,G, JVV)
-C
-C this subroutine computes an off-shell vector current from the three-  
-C point gauge boson coupling.  the vector propagator is given in feynman
-C gauge for a massless vector and in unitary gauge for a massive vector.
-C                                                                       
-C input:                                                                
-C       complex v1(6)          : first  vector                        v1
-C       complex v2(6)          : second vector                        v2
-C       real    g              : coupling constant (see the table below)
-C                                                                       
-C output:                                                               
-C       complex jvv(6)         : vector current            j^mu(v:v1,v2)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 V1(6),V2(6),JVV(6),J12(0:3),
-     &        SV1,SV2,V12
-      REAL*8    P1(0:3),P2(0:3),Q(0:3),G,GS,S
-C
-      REAL*8 RXZERO
-      PARAMETER( RXZERO=0.0D0 )
-C
-      JVV(5) = V1(5)+V2(5)
-      JVV(6) = V1(6)+V2(6)
-C
-      P1(0)=DBLE( V1(5))
-      P1(1)=DBLE( V1(6))
-      P1(2)=DIMAG(V1(6))
-      P1(3)=DIMAG(V1(5))
-      P2(0)=DBLE( V2(5))
-      P2(1)=DBLE( V2(6))
-      P2(2)=DIMAG(V2(6))
-      P2(3)=DIMAG(V2(5))
-      Q(0)=-DBLE( JVV(5))
-      Q(1)=-DBLE( JVV(6))
-      Q(2)=-DIMAG(JVV(6))
-      Q(3)=-DIMAG(JVV(5))
-      S=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-C
-      V12=V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4)
-      SV1= (P2(0)-Q(0))*V1(1) -(P2(1)-Q(1))*V1(2)
-     &    -(P2(2)-Q(2))*V1(3) -(P2(3)-Q(3))*V1(4)
-      SV2=-(P1(0)-Q(0))*V2(1) +(P1(1)-Q(1))*V2(2)
-     &    +(P1(2)-Q(2))*V2(3) +(P1(3)-Q(3))*V2(4)
-      J12(0)=(P1(0)-P2(0))*V12 +SV1*V2(1) +SV2*V1(1)
-      J12(1)=(P1(1)-P2(1))*V12 +SV1*V2(2) +SV2*V1(2)
-      J12(2)=(P1(2)-P2(2))*V12 +SV1*V2(3) +SV2*V1(3)
-      J12(3)=(P1(3)-P2(3))*V12 +SV1*V2(4) +SV2*V1(4)
-C
-      GS=-G/S
-C
-      JVV(1) = GS*J12(0)
-      JVV(2) = GS*J12(1)
-      JVV(3) = GS*J12(2)
-      JVV(4) = GS*J12(3)
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JIOXXX(FI,FO,G,VMASS,VWIDTH , JIO)
-C
-C this subroutine computes an off-shell vector current from an external 
-C fermion pair.  the vector boson propagator is given in feynman gauge  
-C for a massless vector and in unitary gauge for a massive vector.      
-C                                                                       
-C input:                                                                
-C       complex fi(6)          : flow-in  fermion                   |fi>
-C       complex fo(6)          : flow-out fermion                   <fo|
-C       real    g(2)           : coupling constants                  gvf
-C       real    vmass          : mass  of output vector v               
-C       real    vwidth         : width of output vector v               
-C                                                                       
-C output:                                                               
-C       complex jio(6)         : vector current          j^mu(<fo|v|fi>)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FI(6),FO(6),JIO(6),C0,C1,C2,C3,CS,D
-      REAL*8    G(2),Q(0:3),VMASS,VWIDTH,Q2,VM2,DD
-C
-      REAL*8 RXZERO, RXONE
-      PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
-      COMPLEX*16 CXIMAG
-      LOGICAL FIRST
-      SAVE CXIMAG,FIRST
-      DATA FIRST/.TRUE./
-C
-C          Fix compilation with g77
-      IF(FIRST) THEN
-        FIRST=.FALSE.
-        CXIMAG=DCMPLX( RXZERO, RXONE )
-      ENDIF
-C
-      JIO(5) = FO(5)-FI(5)
-      JIO(6) = FO(6)-FI(6)
-C
-      Q(0)=DBLE( JIO(5))
-      Q(1)=DBLE( JIO(6))
-      Q(2)=DIMAG(JIO(6))
-      Q(3)=DIMAG(JIO(5))
-      Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-      VM2=VMASS**2
-C
-      IF (VMASS.NE.RXZERO) THEN
-C
-         D=RXONE/DCMPLX( Q2-VM2 , MAX(SIGN( VMASS*VWIDTH ,Q2),RXZERO) )
-C  for the running width, use below instead of the above d.
-C      d=r_one/dcmplx( q2-vm2 , max( vwidth*q2/vmass ,r_zero) )
-C
-         IF (G(2).NE.RXZERO) THEN
-C
-            C0=  G(1)*( FO(3)*FI(1)+FO(4)*FI(2))
-     &          +G(2)*( FO(1)*FI(3)+FO(2)*FI(4))
-            C1= -G(1)*( FO(3)*FI(2)+FO(4)*FI(1))
-     &          +G(2)*( FO(1)*FI(4)+FO(2)*FI(3))
-            C2=( G(1)*( FO(3)*FI(2)-FO(4)*FI(1)) 
-     &          +G(2)*(-FO(1)*FI(4)+FO(2)*FI(3)))*CXIMAG
-            C3=  G(1)*(-FO(3)*FI(1)+FO(4)*FI(2))
-     &          +G(2)*( FO(1)*FI(3)-FO(2)*FI(4))
-         ELSE
-C
-            D=D*G(1)
-            C0=  FO(3)*FI(1)+FO(4)*FI(2)
-            C1= -FO(3)*FI(2)-FO(4)*FI(1)
-            C2=( FO(3)*FI(2)-FO(4)*FI(1))*CXIMAG
-            C3= -FO(3)*FI(1)+FO(4)*FI(2)
-         END IF
-C
-         CS=(Q(0)*C0-Q(1)*C1-Q(2)*C2-Q(3)*C3)/VM2
-C
-         JIO(1) = (C0-CS*Q(0))*D
-         JIO(2) = (C1-CS*Q(1))*D
-         JIO(3) = (C2-CS*Q(2))*D
-         JIO(4) = (C3-CS*Q(3))*D
-C
-      ELSE
-         DD=RXONE/Q2
-C
-         IF (G(2).NE.RXZERO) THEN
-            JIO(1) = ( G(1)*( FO(3)*FI(1)+FO(4)*FI(2))
-     &                +G(2)*( FO(1)*FI(3)+FO(2)*FI(4)) )*DD
-            JIO(2) = (-G(1)*( FO(3)*FI(2)+FO(4)*FI(1))
-     &                +G(2)*( FO(1)*FI(4)+FO(2)*FI(3)) )*DD
-            JIO(3) = ( G(1)*( FO(3)*FI(2)-FO(4)*FI(1))
-     &                +G(2)*(-FO(1)*FI(4)+FO(2)*FI(3)))
-     $           *DCMPLX(RXZERO,DD)
-            JIO(4) = ( G(1)*(-FO(3)*FI(1)+FO(4)*FI(2))
-     &                +G(2)*( FO(1)*FI(3)-FO(2)*FI(4)) )*DD
-C
-         ELSE
-            DD=DD*G(1)
-C
-            JIO(1) =  ( FO(3)*FI(1)+FO(4)*FI(2))*DD
-            JIO(2) = -( FO(3)*FI(2)+FO(4)*FI(1))*DD
-            JIO(3) =  ( FO(3)*FI(2)-FO(4)*FI(1))*DCMPLX(RXZERO,DD)
-            JIO(4) =  (-FO(3)*FI(1)+FO(4)*FI(2))*DD
-         END IF
-      END IF
-C
-      RETURN
-      END
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JSSXXX(S1,S2,G,VMASS,VWIDTH , JSS)
-C
-C This subroutine computes an off-shell vector current from the vector- 
-C scalar-scalar coupling.  The coupling is absent in the minimal SM in  
-C unitary gauge.  The propagator is given in Feynman gauge for a        
-C massless vector and in unitary gauge for a massive vector.            
-C                                                                       
-C INPUT:                                                                
-C       complex S1(3)          : first  scalar                        S1
-C       complex S2(3)          : second scalar                        S2
-C       real    G              : coupling constant (S1 charge)          
-C       real    VMASS          : mass  of OUTPUT vector V               
-C       real    VWIDTH         : width of OUTPUT vector V               
-C                                                                       
-C Examples of the coupling constant G for SUSY particles are as follows:
-C   -----------------------------------------------------------         
-C   |    S1    | (Q,I3) of S1  ||   V=A   |   V=Z   |   V=W   |         
-C   -----------------------------------------------------------         
-C   | nu~_L    | (  0  , +1/2) ||   ---   |  GZN(1) |  GWF(1) |         
-C   | e~_L     | ( -1  , -1/2) ||  GAL(1) |  GZL(1) |  GWF(1) |         
-C   | u~_L     | (+2/3 , +1/2) ||  GAU(1) |  GZU(1) |  GWF(1) |         
-C   | d~_L     | (-1/3 , -1/2) ||  GAD(1) |  GZD(1) |  GWF(1) |         
-C   -----------------------------------------------------------         
-C   | e~_R-bar | ( +1  ,  0  ) || -GAL(2) | -GZL(2) | -GWF(2) |         
-C   | u~_R-bar | (-2/3 ,  0  ) || -GAU(2) | -GZU(2) | -GWF(2) |         
-C   | d~_R-bar | (+1/3 ,  0  ) || -GAD(2) | -GZD(2) | -GWF(2) |         
-C   -----------------------------------------------------------         
-C where the S1 charge is defined by the flowing-OUT quantum number.     
-C                                                                       
-C OUTPUT:                                                               
-C       complex JSS(6)         : vector current            J^mu(V:S1,S2)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 S1(3),S2(3),JSS(6),DG,ADG
-      REAL*8  PP(0:3),PA(0:3),Q(0:3),G,VMASS,VWIDTH,Q2,VM2,MP2,MA2,M2D
-C
-      JSS(5) = S1(2)+S2(2)
-      JSS(6) = S1(3)+S2(3)
-C
-      Q(0)=DBLE( JSS(5))
-      Q(1)=DBLE( JSS(6))
-      Q(2)=DIMAG(JSS(6))
-      Q(3)=DIMAG(JSS(5))
-      Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-      VM2=VMASS**2
-C
-      IF (VMASS.EQ.0.) GOTO 10
-C
-      DG=G/DCMPLX( Q2-VM2, MAX(SIGN( VMASS*VWIDTH ,Q2),0.D0))
-C  For the running width, use below instead of the above DG.
-C      DG=G/dCMPLX( Q2-VM2 , MAX( VWIDTH*Q2/VMASS ,0.) )
-C
-      ADG=DG*S1(1)*S2(1)
-C
-      PP(0)=DBLE( S1(2))
-      PP(1)=DBLE( S1(3))
-      PP(2)=DIMAG(S1(3))
-      PP(3)=DIMAG(S1(2))
-      PA(0)=DBLE( S2(2))
-      PA(1)=DBLE( S2(3))
-      PA(2)=DIMAG(S2(3))
-      PA(3)=DIMAG(S2(2))
-      MP2=PP(0)**2-(PP(1)**2+PP(2)**2+PP(3)**2)
-      MA2=PA(0)**2-(PA(1)**2+PA(2)**2+PA(3)**2)
-      M2D=MP2-MA2
-C
-      JSS(1) = ADG*( (PP(0)-PA(0)) - Q(0)*M2D/VM2)
-      JSS(2) = ADG*( (PP(1)-PA(1)) - Q(1)*M2D/VM2)
-      JSS(3) = ADG*( (PP(2)-PA(2)) - Q(2)*M2D/VM2)
-      JSS(4) = ADG*( (PP(3)-PA(3)) - Q(3)*M2D/VM2)
-C
-      RETURN
-C
-  10  ADG=G*S1(1)*S2(1)/Q2
-C
-      JSS(1) = ADG*DBLE( S1(2)-S2(2))
-      JSS(2) = ADG*DBLE( S1(3)-S2(3))
-      JSS(3) = ADG*DIMAG(S1(3)-S2(3))
-      JSS(4) = ADG*DIMAG(S1(2)-S2(2))
-C
-      RETURN
-      END
-C
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JTIOXX(FI,FO,G , JIO)
-C
-C this subroutine computes an off-shell vector current from an external 
-C fermion pair.  the vector boson propagator is not included in this
-C routine.
-C                                                                       
-C input:                                                                
-C       complex fi(6)          : flow-in  fermion                   |fi>
-C       complex fo(6)          : flow-out fermion                   <fo|
-C       real    g(2)           : coupling constants                  gvf
-C                                                                       
-C output:                                                               
-C       complex jio(6)         : vector current          j^mu(<fo|v|fi>)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FI(6),FO(6),JIO(6)
-      REAL*8    G(2)
-C
-      REAL*8 RXZERO, RXONE
-      PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
-      COMPLEX*16 CXIMAG
-      LOGICAL FIRST
-      SAVE CXIMAG,FIRST
-      DATA FIRST/.TRUE./
-C
-C          Fix compilation with g77
-      IF(FIRST) THEN
-        FIRST=.FALSE.
-        CXIMAG=DCMPLX( RXZERO, RXONE )
-      ENDIF
-C
-      JIO(5) = FO(5)-FI(5)
-      JIO(6) = FO(6)-FI(6)
-C
-      IF ( G(2) .NE. RXZERO ) THEN
-         JIO(1) = ( G(1)*( FO(3)*FI(1)+FO(4)*FI(2))
-     &             +G(2)*( FO(1)*FI(3)+FO(2)*FI(4)) )
-         JIO(2) = (-G(1)*( FO(3)*FI(2)+FO(4)*FI(1))
-     &             +G(2)*( FO(1)*FI(4)+FO(2)*FI(3)) )
-         JIO(3) = ( G(1)*( FO(3)*FI(2)-FO(4)*FI(1))
-     &             +G(2)*(-FO(1)*FI(4)+FO(2)*FI(3)) )*CXIMAG
-         JIO(4) = ( G(1)*(-FO(3)*FI(1)+FO(4)*FI(2))
-     &             +G(2)*( FO(1)*FI(3)-FO(2)*FI(4)) )
-C
-      ELSE
-         JIO(1) =  ( FO(3)*FI(1)+FO(4)*FI(2))*G(1)
-         JIO(2) = -( FO(3)*FI(2)+FO(4)*FI(1))*G(1)
-         JIO(3) =  ( FO(3)*FI(2)-FO(4)*FI(1))*DCMPLX(RXZERO,G(1))
-         JIO(4) =  (-FO(3)*FI(1)+FO(4)*FI(2))*G(1)
-      END IF
-C
-      RETURN
-      END
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JVSSXX(VC,S1,S2,G,VMASS,VWIDTH , JVSS)
-C
-C This subroutine computes an off-shell vector current from the vector- 
-C vector-scalar-scalar coupling.  The vector propagator is given in     
-C Feynman gauge for a massless vector and in unitary gauge for a massive
-C vector.                                                               
-C                                                                       
-C INPUT:                                                                
-C       complex VC(6)          : input  vector                        V 
-C       complex S1(3)          : first  scalar                        S1
-C       complex S2(3)          : second scalar                        S2
-C       real    G              : coupling constant                 GVVHH
-C       real    VMASS          : mass  of OUTPUT vector V'              
-C       real    VWIDTH         : width of OUTPUT vector V'              
-C                                                                       
-C OUTPUT:                                                               
-C       complex JVSS(6)        : vector current         J^mu(V':V,S1,S2)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 VC(6),S1(3),S2(3),JVSS(6),DG
-      REAL*8    Q(0:3),G,VMASS,VWIDTH,Q2,VK,VM2
-C
-      JVSS(5) = VC(5)+S1(2)+S2(2)
-      JVSS(6) = VC(6)+S1(3)+S2(3)
-C
-      Q(0)=DBLE( JVSS(5))
-      Q(1)=DBLE( JVSS(6))
-      Q(2)=DIMAG(JVSS(6))
-      Q(3)=DIMAG(JVSS(5))
-      Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-      VM2=VMASS**2
-C
-      IF (VMASS.EQ.0.) GOTO 10
-C
-      DG=G*S1(1)*S2(1)/DCMPLX( Q2-VM2,MAX(SIGN( VMASS*VWIDTH,Q2),0.D0))
-C  For the running width, use below instead of the above DG.
-C      DG=G*S1(1)*S2(1)/CMPLX( Q2-VM2 , MAX( VWIDTH*Q2/VMASS ,0.))
-C
-      VK=(Q(0)*VC(1)-Q(1)*VC(2)-Q(2)*VC(3)-Q(3)*VC(4))/VM2
-C
-      JVSS(1) = DG*(VC(1)-VK*Q(0))
-      JVSS(2) = DG*(VC(2)-VK*Q(1))
-      JVSS(3) = DG*(VC(3)-VK*Q(2))
-      JVSS(4) = DG*(VC(4)-VK*Q(3))
-C
-      RETURN
-C
-  10  DG= G*S1(1)*S2(1)/Q2
-C
-      JVSS(1) = DG*VC(1)
-      JVSS(2) = DG*VC(2)
-      JVSS(3) = DG*VC(3)
-      JVSS(4) = DG*VC(4)
-C
-      RETURN
-      END
-C
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JVSXXX(VC,SC,G,VMASS,VWIDTH , JVS)
-C
-C this subroutine computes an off-shell vector current from the vector- 
-C vector-scalar coupling.  the vector propagator is given in feynman    
-C gauge for a massless vector and in unitary gauge for a massive vector.
-C                                                                       
-C input:                                                                
-C       complex vc(6)          : input vector                          v
-C       complex sc(3)          : input scalar                          s
-C       real    g              : coupling constant                  gvvh
-C       real    vmass          : mass  of output vector v'              
-C       real    vwidth         : width of output vector v'              
-C                                                                       
-C output:                                                               
-C       complex jvs(6)         : vector current             j^mu(v':v,s)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 VC(6),SC(3),JVS(6),DG,VK
-      REAL*8    Q(0:3),VMASS,VWIDTH,Q2,VM2,G
-C
-      JVS(5) = VC(5)+SC(2)
-      JVS(6) = VC(6)+SC(3)
-C
-      Q(0)=DBLE( JVS(5))
-      Q(1)=DBLE( JVS(6))
-      Q(2)=DIMAG(JVS(6))
-      Q(3)=DIMAG(JVS(5))
-      Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-      VM2=VMASS**2
-C
-      IF (VMASS.EQ.0.) GOTO 10
-C
-      DG=G*SC(1)/DCMPLX( Q2-VM2 , MAX(DSIGN( VMASS*VWIDTH ,Q2),0.D0) )
-C  for the running width, use below instead of the above dg.
-C      dg=g*sc(1)/dcmplx( q2-vm2 , max( vwidth*q2/vmass ,0.) )
-C
-      VK=(-Q(0)*VC(1)+Q(1)*VC(2)+Q(2)*VC(3)+Q(3)*VC(4))/VM2
-C
-      JVS(1) = DG*(Q(0)*VK+VC(1))
-      JVS(2) = DG*(Q(1)*VK+VC(2))
-      JVS(3) = DG*(Q(2)*VK+VC(3))
-      JVS(4) = DG*(Q(3)*VK+VC(4))
-C
-      RETURN
-C
-  10  DG=G*SC(1)/Q2
-C
-      JVS(1) = DG*VC(1)
-      JVS(2) = DG*VC(2)
-      JVS(3) = DG*VC(3)
-      JVS(4) = DG*VC(4)
-C
-      RETURN
-      END
-
-
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JVVXXX(V1,V2,G,VMASS,VWIDTH , JVV)
-C
-C this subroutine computes an off-shell vector current from the three-  
-C point gauge boson coupling.  the vector propagator is given in feynman
-C gauge for a massless vector and in unitary gauge for a massive vector.
-C                                                                       
-C input:                                                                
-C       complex v1(6)          : first  vector                        v1
-C       complex v2(6)          : second vector                        v2
-C       real    g              : coupling constant (see the table below)
-C       real    vmass          : mass  of output vector v               
-C       real    vwidth         : width of output vector v               
-C                                                                       
-C the possible sets of the inputs are as follows:                       
-C    ------------------------------------------------------------------ 
-C    |   v1   |   v2   |  jvv   |      g       |   vmass  |  vwidth   | 
-C    ------------------------------------------------------------------ 
-C    |   w-   |   w+   |  a/z   |  gwwa/gwwz   | 0./zmass | 0./zwidth | 
-C    | w3/a/z |   w-   |  w+    | gw/gwwa/gwwz |   wmass  |  wwidth   | 
-C    |   w+   | w3/a/z |  w-    | gw/gwwa/gwwz |   wmass  |  wwidth   | 
-C    ------------------------------------------------------------------ 
-C where all the bosons are defined by the flowing-out quantum number.   
-C                                                                       
-C output:                                                               
-C       complex jvv(6)         : vector current            j^mu(v:v1,v2)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 V1(6),V2(6),JVV(6),J12(0:3),JS,DG,
-     &        SV1,SV2,S11,S12,S21,S22,V12
-      REAL*8    P1(0:3),P2(0:3),Q(0:3),G,VMASS,VWIDTH,GS,S,VM2,M1,M2
-C
-      REAL*8 RXZERO
-      PARAMETER( RXZERO=0.0D0 )
-C
-      JVV(5) = V1(5)+V2(5)
-      JVV(6) = V1(6)+V2(6)
-C
-      P1(0)=DBLE( V1(5))
-      P1(1)=DBLE( V1(6))
-      P1(2)=DIMAG(V1(6))
-      P1(3)=DIMAG(V1(5))
-      P2(0)=DBLE( V2(5))
-      P2(1)=DBLE( V2(6))
-      P2(2)=DIMAG(V2(6))
-      P2(3)=DIMAG(V2(5))
-      Q(0)=-DBLE( JVV(5))
-      Q(1)=-DBLE( JVV(6))
-      Q(2)=-DIMAG(JVV(6))
-      Q(3)=-DIMAG(JVV(5))
-      S=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
-C
-      V12=V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4)
-      SV1= (P2(0)-Q(0))*V1(1) -(P2(1)-Q(1))*V1(2)
-     &    -(P2(2)-Q(2))*V1(3) -(P2(3)-Q(3))*V1(4)
-      SV2=-(P1(0)-Q(0))*V2(1) +(P1(1)-Q(1))*V2(2)
-     &    +(P1(2)-Q(2))*V2(3) +(P1(3)-Q(3))*V2(4)
-      J12(0)=(P1(0)-P2(0))*V12 +SV1*V2(1) +SV2*V1(1)
-      J12(1)=(P1(1)-P2(1))*V12 +SV1*V2(2) +SV2*V1(2)
-      J12(2)=(P1(2)-P2(2))*V12 +SV1*V2(3) +SV2*V1(3)
-      J12(3)=(P1(3)-P2(3))*V12 +SV1*V2(4) +SV2*V1(4)
-C
-      IF ( VMASS .NE. RXZERO ) THEN
-         VM2=VMASS**2
-         M1=P1(0)**2-(P1(1)**2+P1(2)**2+P1(3)**2)
-         M2=P2(0)**2-(P2(1)**2+P2(2)**2+P2(3)**2)
-         S11=P1(0)*V1(1)-P1(1)*V1(2)-P1(2)*V1(3)-P1(3)*V1(4)
-         S12=P1(0)*V2(1)-P1(1)*V2(2)-P1(2)*V2(3)-P1(3)*V2(4)
-         S21=P2(0)*V1(1)-P2(1)*V1(2)-P2(2)*V1(3)-P2(3)*V1(4)
-         S22=P2(0)*V2(1)-P2(1)*V2(2)-P2(2)*V2(3)-P2(3)*V2(4)
-         JS=(V12*(-M1+M2) +S11*S12 -S21*S22)/VM2
-C
-         DG=-G/DCMPLX( S-VM2 , MAX(SIGN( VMASS*VWIDTH ,S),RXZERO) )
-C
-C  for the running width, use below instead of the above dg.
-C         dg=-g/dcmplx( s-vm2 , max( vwidth*s/vmass ,r_zero) )
-C
-         JVV(1) = DG*(J12(0)-Q(0)*JS)
-         JVV(2) = DG*(J12(1)-Q(1)*JS)
-         JVV(3) = DG*(J12(2)-Q(2)*JS)
-         JVV(4) = DG*(J12(3)-Q(3)*JS)
-C
-      ELSE
-         GS=-G/S
-C
-         JVV(1) = GS*J12(0)
-         JVV(2) = GS*J12(1)
-         JVV(3) = GS*J12(2)
-         JVV(4) = GS*J12(3)
-      END IF
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JW3WXX(W1,W2,W3,G1,G2,WMASS,WWIDTH,VMASS,VWIDTH , JW3W)
-C
-C this subroutine computes an off-shell w+, w-, w3, z or photon current 
-C from the four-point gauge boson coupling, including the contributions 
-C of w exchange diagrams.  the vector propagator is given in feynman    
-C gauge for a photon and in unitary gauge for w and z bosons.  if one   
-C sets wmass=0.0, then the ggg-->g current is given (see sect 2.9.1 of 
-C the manual).                                                          
-C                                                                       
-C input:                                                                
-C       complex w1(6)          : first  vector                        w1
-C       complex w2(6)          : second vector                        w2
-C       complex w3(6)          : third  vector                        w3
-C       real    g1             : first  coupling constant               
-C       real    g2             : second coupling constant               
-C                                                  (see the table below)
-C       real    wmass          : mass  of internal w                    
-C       real    wwidth         : width of internal w                    
-C       real    vmass          : mass  of output w'                     
-C       real    vwidth         : width of output w'                     
-C                                                                       
-C the possible sets of the inputs are as follows:                       
-C   ------------------------------------------------------------------- 
-C   |  w1  |  w2  |  w3  | g1 | g2 |wmass|wwidth|vmass|vwidth || jw3w | 
-C   ------------------------------------------------------------------- 
-C   |  w-  |  w3  |  w+  | gw |gwwz|wmass|wwidth|zmass|zwidth ||  z   | 
-C   |  w-  |  w3  |  w+  | gw |gwwa|wmass|wwidth|  0. |  0.   ||  a   | 
-C   |  w-  |  z   |  w+  |gwwz|gwwz|wmass|wwidth|zmass|zwidth ||  z   | 
-C   |  w-  |  z   |  w+  |gwwz|gwwa|wmass|wwidth|  0. |  0.   ||  a   | 
-C   |  w-  |  a   |  w+  |gwwa|gwwz|wmass|wwidth|zmass|zwidth ||  z   | 
-C   |  w-  |  a   |  w+  |gwwa|gwwa|wmass|wwidth|  0. |  0.   ||  a   | 
-C   ------------------------------------------------------------------- 
-C   |  w3  |  w-  |  w3  | gw | gw |wmass|wwidth|wmass|wwidth ||  w+  | 
-C   |  w3  |  w+  |  w3  | gw | gw |wmass|wwidth|wmass|wwidth ||  w-  | 
-C   |  w3  |  w-  |  z   | gw |gwwz|wmass|wwidth|wmass|wwidth ||  w+  | 
-C   |  w3  |  w+  |  z   | gw |gwwz|wmass|wwidth|wmass|wwidth ||  w-  | 
-C   |  w3  |  w-  |  a   | gw |gwwa|wmass|wwidth|wmass|wwidth ||  w+  | 
-C   |  w3  |  w+  |  a   | gw |gwwa|wmass|wwidth|wmass|wwidth ||  w-  | 
-C   |  z   |  w-  |  z   |gwwz|gwwz|wmass|wwidth|wmass|wwidth ||  w+  | 
-C   |  z   |  w+  |  z   |gwwz|gwwz|wmass|wwidth|wmass|wwidth ||  w-  | 
-C   |  z   |  w-  |  a   |gwwz|gwwa|wmass|wwidth|wmass|wwidth ||  w+  | 
-C   |  z   |  w+  |  a   |gwwz|gwwa|wmass|wwidth|wmass|wwidth ||  w-  | 
-C   |  a   |  w-  |  a   |gwwa|gwwa|wmass|wwidth|wmass|wwidth ||  w+  | 
-C   |  a   |  w+  |  a   |gwwa|gwwa|wmass|wwidth|wmass|wwidth ||  w-  | 
-C   ------------------------------------------------------------------- 
-C where all the bosons are defined by the flowing-out quantum number.   
-C                                                                       
-C output:                                                               
-C       complex jw3w(6)        : w current             j^mu(w':w1,w2,w3)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16  W1(6),W2(6),W3(6),JW3W(6)
-      COMPLEX*16 DW1(0:3),DW2(0:3),DW3(0:3),
-     &           JJ(0:3),J4(0:3),
-     &           DV,W12,W32,W13,
-     &           JQ
-      REAL*8     G1,G2,WMASS,WWIDTH,VMASS,VWIDTH
-      REAL*8     P1(0:3),P2(0:3),P3(0:3),Q(0:3),
-     &           DG2,DMV,DWV,MV2,Q2
-C
-      REAL*8 RXZERO
-      PARAMETER( RXZERO=0.0D0 )
-C
-      JW3W(5) = W1(5)+W2(5)+W3(5)
-      JW3W(6) = W1(6)+W2(6)+W3(6)
-C
-      DW1(0)=DCMPLX(W1(1))
-      DW1(1)=DCMPLX(W1(2))
-      DW1(2)=DCMPLX(W1(3))
-      DW1(3)=DCMPLX(W1(4))
-      DW2(0)=DCMPLX(W2(1))
-      DW2(1)=DCMPLX(W2(2))
-      DW2(2)=DCMPLX(W2(3))
-      DW2(3)=DCMPLX(W2(4))
-      DW3(0)=DCMPLX(W3(1))
-      DW3(1)=DCMPLX(W3(2))
-      DW3(2)=DCMPLX(W3(3))
-      DW3(3)=DCMPLX(W3(4))
-      P1(0)=DBLE(      W1(5))
-      P1(1)=DBLE(      W1(6))
-      P1(2)=DBLE(DIMAG(W1(6)))
-      P1(3)=DBLE(DIMAG(W1(5)))
-      P2(0)=DBLE(      W2(5))
-      P2(1)=DBLE(      W2(6))
-      P2(2)=DBLE(DIMAG(W2(6)))
-      P2(3)=DBLE(DIMAG(W2(5)))
-      P3(0)=DBLE(      W3(5))
-      P3(1)=DBLE(      W3(6))
-      P3(2)=DBLE(DIMAG(W3(6)))
-      P3(3)=DBLE(DIMAG(W3(5)))
-      Q(0)=-(P1(0)+P2(0)+P3(0))
-      Q(1)=-(P1(1)+P2(1)+P3(1))
-      Q(2)=-(P1(2)+P2(2)+P3(2))
-      Q(3)=-(P1(3)+P2(3)+P3(3))
-
-
-      Q2 =Q(0)**2 -(Q(1)**2 +Q(2)**2 +Q(3)**2)
-      DG2=DBLE(G1)*DBLE(G2)
-      DMV=DBLE(VMASS)
-      DWV=DBLE(VWIDTH)
-      MV2=DMV**2
-      IF (VMASS.EQ. RXZERO) THEN
-      DV = 1.0D0/DCMPLX( Q2 )
-      ELSE
-      DV = 1.0D0/DCMPLX( Q2 -MV2 , DMAX1(DSIGN(DMV*DWV,Q2 ),0.D0) )
-      ENDIF
-C  for the running width, use below instead of the above dv.
-C      dv = 1.0d0/dcmplx( q2 -mv2 , dmax1(dwv*q2/dmv,0.d0) )
-C
-      W12=DW1(0)*DW2(0)-DW1(1)*DW2(1)-DW1(2)*DW2(2)-DW1(3)*DW2(3)
-      W32=DW3(0)*DW2(0)-DW3(1)*DW2(1)-DW3(2)*DW2(2)-DW3(3)*DW2(3)
-C
-      IF ( WMASS .NE. RXZERO ) THEN
-         W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
-C
-         J4(0)=DG2*( DW1(0)*W32 + DW3(0)*W12 - 2.D0*DW2(0)*W13 )
-         J4(1)=DG2*( DW1(1)*W32 + DW3(1)*W12 - 2.D0*DW2(1)*W13 )
-         J4(2)=DG2*( DW1(2)*W32 + DW3(2)*W12 - 2.D0*DW2(2)*W13 )
-         J4(3)=DG2*( DW1(3)*W32 + DW3(3)*W12 - 2.D0*DW2(3)*W13 )
-C
-         JJ(0)=J4(0)
-         JJ(1)=J4(1)
-         JJ(2)=J4(2)
-         JJ(3)=J4(3)
-
-      ELSE
-C
-         W12=DW1(0)*DW2(0)-DW1(1)*DW2(1)-DW1(2)*DW2(2)-DW1(3)*DW2(3)
-         W32=DW3(0)*DW2(0)-DW3(1)*DW2(1)-DW3(2)*DW2(2)-DW3(3)*DW2(3)
-         W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
-C
-         J4(0)=DG2*( DW1(0)*W32 - DW2(0)*W13 )
-         J4(1)=DG2*( DW1(1)*W32 - DW2(1)*W13 )
-         J4(2)=DG2*( DW1(2)*W32 - DW2(2)*W13 )
-         J4(3)=DG2*( DW1(3)*W32 - DW2(3)*W13 )
-C
-         JJ(0)=J4(0)
-         JJ(1)=J4(1)
-         JJ(2)=J4(2)
-         JJ(3)=J4(3)
-
-      END IF
-C
-      IF ( VMASS .NE. RXZERO ) THEN
-C
-         JQ=(JJ(0)*Q(0)-JJ(1)*Q(1)-JJ(2)*Q(2)-JJ(3)*Q(3))/MV2
-C
-         JW3W(1) = DCMPLX( (JJ(0)-JQ*Q(0))*DV )
-         JW3W(2) = DCMPLX( (JJ(1)-JQ*Q(1))*DV )
-         JW3W(3) = DCMPLX( (JJ(2)-JQ*Q(2))*DV )
-         JW3W(4) = DCMPLX( (JJ(3)-JQ*Q(3))*DV )
-C
-      ELSE
-C
-         JW3W(1) = DCMPLX( JJ(0)*DV )
-         JW3W(2) = DCMPLX( JJ(1)*DV )
-         JW3W(3) = DCMPLX( JJ(2)*DV )
-         JW3W(4) = DCMPLX( JJ(3)*DV )
-      END IF
-C
-      RETURN
-      END
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE JWWWXX(W1,W2,W3,GWWA,GWWZ,ZMASS,ZWIDTH,WMASS,WWIDTH ,
-     &                  JWWW)
-C
-C this subroutine computes an off-shell w+/w- current from the four-    
-C point gauge boson coupling, including the contributions of photon and 
-C z exchanges.  the vector propagators for the output w and the internal
-C z bosons are given in unitary gauge, and that of the internal photon  
-C is given in feynman gauge.                                            
-C                                                                       
-C input:                                                                
-C       complex w1(6)          : first  vector                        w1
-C       complex w2(6)          : second vector                        w2
-C       complex w3(6)          : third  vector                        w3
-C       real    gwwa           : coupling constant of w and a       gwwa
-C       real    gwwz           : coupling constant of w and z       gwwz
-C       real    zmass          : mass  of internal z                    
-C       real    zwidth         : width of internal z                    
-C       real    wmass          : mass  of output w                      
-C       real    wwidth         : width of output w                      
-C                                                                       
-C the possible sets of the inputs are as follows:                       
-C   ------------------------------------------------------------------- 
-C   |  w1  |  w2  |  w3  |gwwa|gwwz|zmass|zwidth|wmass|wwidth || jwww | 
-C   ------------------------------------------------------------------- 
-C   |  w-  |  w+  |  w-  |gwwa|gwwz|zmass|zwidth|wmass|wwidth ||  w+  | 
-C   |  w+  |  w-  |  w+  |gwwa|gwwz|zmass|zwidth|wmass|wwidth ||  w-  | 
-C   ------------------------------------------------------------------- 
-C where all the bosons are defined by the flowing-out quantum number.   
-C                                                                       
-C output:                                                               
-C       complex jwww(6)        : w current             j^mu(w':w1,w2,w3)
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16  W1(6),W2(6),W3(6),JWWW(6)
-      COMPLEX*16 DW1(0:3),DW2(0:3),DW3(0:3),
-     &           JJ(0:3),JS(0:3),JT(0:3),J4(0:3),
-     &           JT12(0:3),JT32(0:3),J12(0:3),J32(0:3),
-     &           DZS,DZT,DW,W12,W32,W13,P1W2,P2W1,P3W2,P2W3,
-     &           JK12,JK32,JSW3,JTW1,P3JS,KSW3,P1JT,KTW1,JQ
-      REAL*8     GWWA,GWWZ,ZMASS,ZWIDTH,WMASS,WWIDTH
-      REAL*8     P1(0:3),P2(0:3),P3(0:3),Q(0:3),KS(0:3),KT(0:3),
-     &           DGWWA2,DGWWZ2,DGW2,DMZ,DWZ,DMW,DWW,MZ2,MW2,Q2,KS2,KT2,
-     &           DAS,DAT
-C
-      JWWW(5) = W1(5)+W2(5)+W3(5)
-      JWWW(6) = W1(6)+W2(6)+W3(6)
-C
-      DW1(0)=DCMPLX(W1(1))
-      DW1(1)=DCMPLX(W1(2))
-      DW1(2)=DCMPLX(W1(3))
-      DW1(3)=DCMPLX(W1(4))
-      DW2(0)=DCMPLX(W2(1))
-      DW2(1)=DCMPLX(W2(2))
-      DW2(2)=DCMPLX(W2(3))
-      DW2(3)=DCMPLX(W2(4))
-      DW3(0)=DCMPLX(W3(1))
-      DW3(1)=DCMPLX(W3(2))
-      DW3(2)=DCMPLX(W3(3))
-      DW3(3)=DCMPLX(W3(4))
-      P1(0)=DBLE(      W1(5))
-      P1(1)=DBLE(      W1(6))
-      P1(2)=DBLE(DIMAG(W1(6)))
-      P1(3)=DBLE(DIMAG(W1(5)))
-      P2(0)=DBLE(      W2(5))
-      P2(1)=DBLE(      W2(6))
-      P2(2)=DBLE(DIMAG(W2(6)))
-      P2(3)=DBLE(DIMAG(W2(5)))
-      P3(0)=DBLE(      W3(5))
-      P3(1)=DBLE(      W3(6))
-      P3(2)=DBLE(DIMAG(W3(6)))
-      P3(3)=DBLE(DIMAG(W3(5)))
-      Q(0)=-(P1(0)+P2(0)+P3(0))
-      Q(1)=-(P1(1)+P2(1)+P3(1))
-      Q(2)=-(P1(2)+P2(2)+P3(2))
-      Q(3)=-(P1(3)+P2(3)+P3(3))
-      KS(0)=P1(0)+P2(0)
-      KS(1)=P1(1)+P2(1)
-      KS(2)=P1(2)+P2(2)
-      KS(3)=P1(3)+P2(3)
-      KT(0)=P2(0)+P3(0)
-      KT(1)=P2(1)+P3(1)
-      KT(2)=P2(2)+P3(2)
-      KT(3)=P2(3)+P3(3)
-      Q2 =Q(0)**2 -(Q(1)**2 +Q(2)**2 +Q(3)**2)
-      KS2=KS(0)**2-(KS(1)**2+KS(2)**2+KS(3)**2)
-      KT2=KT(0)**2-(KT(1)**2+KT(2)**2+KT(3)**2)
-      DGWWA2=DBLE(GWWA)**2
-      DGWWZ2=DBLE(GWWZ)**2
-      DGW2  =DGWWA2+DGWWZ2
-      DMZ=DBLE(ZMASS)
-      DWZ=DBLE(ZWIDTH)
-      DMW=DBLE(WMASS)
-      DWW=DBLE(WWIDTH)
-      MZ2=DMZ**2
-      MW2=DMW**2
-C
-      DAS=-DGWWA2/KS2
-      DAT=-DGWWA2/KT2
-      DZS=-DGWWZ2/DCMPLX( KS2-MZ2 , DMAX1(DSIGN(DMZ*DWZ,KS2),0.D0) )
-      DZT=-DGWWZ2/DCMPLX( KT2-MZ2 , DMAX1(DSIGN(DMZ*DWZ,KT2),0.D0) )
-      DW =-1.0D0/DCMPLX( Q2 -MW2 , DMAX1(DSIGN(DMW*DWW,Q2 ),0.D0) )
-C  for the running width, use below instead of the above dw.
-C      dw =-1.0d0/dcmplx( q2 -mw2 , dmax1(dww*q2/dmw,0.d0) )
-C
-      W12=DW1(0)*DW2(0)-DW1(1)*DW2(1)-DW1(2)*DW2(2)-DW1(3)*DW2(3)
-      W32=DW3(0)*DW2(0)-DW3(1)*DW2(1)-DW3(2)*DW2(2)-DW3(3)*DW2(3)
-C
-      P1W2= (P1(0)+KS(0))*DW2(0)-(P1(1)+KS(1))*DW2(1)
-     &     -(P1(2)+KS(2))*DW2(2)-(P1(3)+KS(3))*DW2(3)
-      P2W1= (P2(0)+KS(0))*DW1(0)-(P2(1)+KS(1))*DW1(1)
-     &     -(P2(2)+KS(2))*DW1(2)-(P2(3)+KS(3))*DW1(3)
-      P3W2= (P3(0)+KT(0))*DW2(0)-(P3(1)+KT(1))*DW2(1)
-     &     -(P3(2)+KT(2))*DW2(2)-(P3(3)+KT(3))*DW2(3)
-      P2W3= (P2(0)+KT(0))*DW3(0)-(P2(1)+KT(1))*DW3(1)
-     &     -(P2(2)+KT(2))*DW3(2)-(P2(3)+KT(3))*DW3(3)
-C
-      JT12(0)= (P1(0)-P2(0))*W12 + P2W1*DW2(0) - P1W2*DW1(0)
-      JT12(1)= (P1(1)-P2(1))*W12 + P2W1*DW2(1) - P1W2*DW1(1)
-      JT12(2)= (P1(2)-P2(2))*W12 + P2W1*DW2(2) - P1W2*DW1(2)
-      JT12(3)= (P1(3)-P2(3))*W12 + P2W1*DW2(3) - P1W2*DW1(3)
-      JT32(0)= (P3(0)-P2(0))*W32 + P2W3*DW2(0) - P3W2*DW3(0)
-      JT32(1)= (P3(1)-P2(1))*W32 + P2W3*DW2(1) - P3W2*DW3(1)
-      JT32(2)= (P3(2)-P2(2))*W32 + P2W3*DW2(2) - P3W2*DW3(2)
-      JT32(3)= (P3(3)-P2(3))*W32 + P2W3*DW2(3) - P3W2*DW3(3)
-C
-      JK12=(JT12(0)*KS(0)-JT12(1)*KS(1)-JT12(2)*KS(2)-JT12(3)*KS(3))/MZ2
-      JK32=(JT32(0)*KT(0)-JT32(1)*KT(1)-JT32(2)*KT(2)-JT32(3)*KT(3))/MZ2
-C
-      J12(0)=JT12(0)*(DAS+DZS)-KS(0)*JK12*DZS
-      J12(1)=JT12(1)*(DAS+DZS)-KS(1)*JK12*DZS
-      J12(2)=JT12(2)*(DAS+DZS)-KS(2)*JK12*DZS
-      J12(3)=JT12(3)*(DAS+DZS)-KS(3)*JK12*DZS
-      J32(0)=JT32(0)*(DAT+DZT)-KT(0)*JK32*DZT
-      J32(1)=JT32(1)*(DAT+DZT)-KT(1)*JK32*DZT
-      J32(2)=JT32(2)*(DAT+DZT)-KT(2)*JK32*DZT
-      J32(3)=JT32(3)*(DAT+DZT)-KT(3)*JK32*DZT
-C
-      JSW3=J12(0)*DW3(0)-J12(1)*DW3(1)-J12(2)*DW3(2)-J12(3)*DW3(3)
-      JTW1=J32(0)*DW1(0)-J32(1)*DW1(1)-J32(2)*DW1(2)-J32(3)*DW1(3)
-C
-      P3JS= (P3(0)-Q(0))*J12(0)-(P3(1)-Q(1))*J12(1)
-     &     -(P3(2)-Q(2))*J12(2)-(P3(3)-Q(3))*J12(3)
-      KSW3= (KS(0)-Q(0))*DW3(0)-(KS(1)-Q(1))*DW3(1)
-     &     -(KS(2)-Q(2))*DW3(2)-(KS(3)-Q(3))*DW3(3)
-      P1JT= (P1(0)-Q(0))*J32(0)-(P1(1)-Q(1))*J32(1)
-     &     -(P1(2)-Q(2))*J32(2)-(P1(3)-Q(3))*J32(3)
-      KTW1= (KT(0)-Q(0))*DW1(0)-(KT(1)-Q(1))*DW1(1)
-     &     -(KT(2)-Q(2))*DW1(2)-(KT(3)-Q(3))*DW1(3)
-C
-      JS(0)= (KS(0)-P3(0))*JSW3 + P3JS*DW3(0) - KSW3*J12(0)
-      JS(1)= (KS(1)-P3(1))*JSW3 + P3JS*DW3(1) - KSW3*J12(1)
-      JS(2)= (KS(2)-P3(2))*JSW3 + P3JS*DW3(2) - KSW3*J12(2)
-      JS(3)= (KS(3)-P3(3))*JSW3 + P3JS*DW3(3) - KSW3*J12(3)
-      JT(0)= (KT(0)-P1(0))*JTW1 + P1JT*DW1(0) - KTW1*J32(0)
-      JT(1)= (KT(1)-P1(1))*JTW1 + P1JT*DW1(1) - KTW1*J32(1)
-      JT(2)= (KT(2)-P1(2))*JTW1 + P1JT*DW1(2) - KTW1*J32(2)
-      JT(3)= (KT(3)-P1(3))*JTW1 + P1JT*DW1(3) - KTW1*J32(3)
-C
-      W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
-C
-      J4(0)=DGW2*( DW1(0)*W32 + DW3(0)*W12 - 2.D0*DW2(0)*W13 )
-      J4(1)=DGW2*( DW1(1)*W32 + DW3(1)*W12 - 2.D0*DW2(1)*W13 )
-      J4(2)=DGW2*( DW1(2)*W32 + DW3(2)*W12 - 2.D0*DW2(2)*W13 )
-      J4(3)=DGW2*( DW1(3)*W32 + DW3(3)*W12 - 2.D0*DW2(3)*W13 )
-C
-C      jj(0)=js(0)+jt(0)+j4(0)
-C      jj(1)=js(1)+jt(1)+j4(1)
-C      jj(2)=js(2)+jt(2)+j4(2)
-C      jj(3)=js(3)+jt(3)+j4(3)
-
-      JJ(0)=J4(0)
-      JJ(1)=J4(1)
-      JJ(2)=J4(2)
-      JJ(3)=J4(3)
-C
-      JQ=(JJ(0)*Q(0)-JJ(1)*Q(1)-JJ(2)*Q(2)-JJ(3)*Q(3))/MW2
-C
-
-      JWWW(1) = DCMPLX( (JJ(0)-JQ*Q(0))*DW )
-      JWWW(2) = DCMPLX( (JJ(1)-JQ*Q(1))*DW )
-      JWWW(3) = DCMPLX( (JJ(2)-JQ*Q(2))*DW )
-      JWWW(4) = DCMPLX( (JJ(3)-JQ*Q(3))*DW )
-C
-      RETURN
-      END
-
-C
-C ----------------------------------------------------------------------
-C
-      SUBROUTINE MOM2CX(ESUM,MASS1,MASS2,COSTH1,PHI1 , P1,P2)
-C
-C This subroutine sets up two four-momenta in the two particle rest     
-C frame.                                                                
-C                                                                       
-C INPUT:                                                                
-C       real    ESUM           : energy sum of particle 1 and 2         
-C       real    MASS1          : mass            of particle 1          
-C       real    MASS2          : mass            of particle 2          
-C       real    COSTH1         : cos(theta)      of particle 1          
-C       real    PHI1           : azimuthal angle of particle 1          
-C                                                                       
-C OUTPUT:                                                               
-C       real    P1(0:3)        : four-momentum of particle 1            
-C       real    P2(0:3)        : four-momentum of particle 2            
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      REAL*8    P1(0:3),P2(0:3),
-     &        ESUM,MASS1,MASS2,COSTH1,PHI1,MD2,ED,PP,SINTH1
-C
-      MD2=(MASS1-MASS2)*(MASS1+MASS2)
-      ED=MD2/ESUM
-      IF (MASS1*MASS2.EQ.0.) THEN
-      PP=(ESUM-ABS(ED))*0.5D0
-C
-      ELSE
-      PP=SQRT((MD2/ESUM)**2-2.0D0*(MASS1**2+MASS2**2)+ESUM**2)*0.5D0
-      ENDIF
-      SINTH1=SQRT((1.0D0-COSTH1)*(1.0D0+COSTH1))
-C
-      P1(0) = MAX((ESUM+ED)*0.5D0,0.D0)
-      P1(1) = PP*SINTH1*COS(PHI1)
-      P1(2) = PP*SINTH1*SIN(PHI1)
-      P1(3) = PP*COSTH1
-C
-      P2(0) = MAX((ESUM-ED)*0.5D0,0.D0)
-      P2(1) = -P1(1)
-      P2(2) = -P1(2)
-      P2(3) = -P1(3)
-C
-      RETURN
-      END
-C **********************************************************************
-C
-      SUBROUTINE MOMNTX(ENERGY,MASS,COSTH,PHI , P)
-C
-C This subroutine sets up a four-momentum from the four inputs.         
-C                                                                       
-C INPUT:                                                                
-C       real    ENERGY         : energy                                 
-C       real    MASS           : mass                                   
-C       real    COSTH          : cos(theta)                             
-C       real    PHI            : azimuthal angle                        
-C                                                                       
-C OUTPUT:                                                               
-C       real    P(0:3)         : four-momentum                          
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      REAL*8    P(0:3),ENERGY,MASS,COSTH,PHI,PP,SINTH
-C
-      P(0) = ENERGY
-      IF (ENERGY.EQ.MASS) THEN
-         P(1) = 0.
-         P(2) = 0.
-         P(3) = 0.
-      ELSE
-         PP=SQRT((ENERGY-MASS)*(ENERGY+MASS))
-         SINTH=SQRT((1.-COSTH)*(1.+COSTH))
-         P(3) = PP*COSTH
-         IF (PHI.EQ.0.) THEN
-            P(1) = PP*SINTH
-            P(2) = 0.
-         ELSE
-            P(1) = PP*SINTH*COS(PHI)
-            P(2) = PP*SINTH*SIN(PHI)
-         ENDIF
-      ENDIF
-      RETURN
-      END
-C
-C
-C
-C      Subroutine returns the desired fermion or
-C      anti-fermion anti-spinor. ie., <f|
-C      A replacement for the HELAS routine OXXXXX
-C
-C      Adam Duff,  1992 August 31
-C      <duff@phenom.physics.wisc.edu>
-C
-      SUBROUTINE OXXXXX(P,FMASS,NHEL,NSF,FO)
-C
-C          P            IN: FOUR VECTOR MOMENTUM
-C          FMASS        IN: FERMION MASS
-C          NHEL         IN: ANTI-SPINOR HELICITY, -1 OR 1
-C          NSF          IN: -1=ANTIFERMION, 1=FERMION
-C          FO           OUT: FERMION WAVEFUNCTION
-C
-C declare input/output variables
-C
-#if defined(CERNLIB_IMPNONE)
-      IMPLICIT NONE
-#endif
-      COMPLEX*16 FO(6)
-      INTEGER*4 NHEL, NSF
-      REAL*8 P(0:3), FMASS
-C
-C declare local variables
-C
-      REAL*8 RXZERO, RXONE, RXTWO
-      PARAMETER( RXZERO=0.0D0, RXONE=1.0D0, RXTWO=2.0D0 )
-      REAL*8 PLAT, PABS, OMEGAP, OMEGAM, RS2PA, SPAZ
-      COMPLEX*16 CXZERO
-      LOGICAL FIRST
-      SAVE CXZERO,FIRST
-      DATA FIRST/.TRUE./
-C
-C          Fix compilation with g77
-      IF(FIRST) THEN
-        FIRST=.FALSE.
-        CXZERO=DCMPLX( RXZERO, RXZERO )
-      ENDIF
-C
-C define kinematic parameters
-C
-      FO(5) = DCMPLX( P(0), P(3) ) * NSF
-      FO(6) = DCMPLX( P(1), P(2) ) * NSF
-      PLAT = SQRT( P(1)**2 + P(2)**2 )
-      PABS = SQRT( P(1)**2 + P(2)**2 + P(3)**2 )
-      OMEGAP = SQRT( P(0) + PABS )
-C
-C do massive fermion case
-C
-      IF ( FMASS .NE. RXZERO ) THEN
-         OMEGAM = FMASS / OMEGAP
-         IF ( NSF .EQ. 1 ) THEN
-            IF ( NHEL .EQ. 1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FO(1) = DCMPLX( OMEGAP, RXZERO )
-                     FO(2) = CXZERO
-                     FO(3) = DCMPLX( OMEGAM, RXZERO )
-                     FO(4) = CXZERO
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS + P(3) )
-                     FO(1) = OMEGAP * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FO(2) = OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( P(1), -P(2) )
-                     FO(3) = OMEGAM * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FO(4) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( P(1), -P(2) )
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FO(1) = CXZERO
-                     FO(2) = DCMPLX( OMEGAP, RXZERO )
-                     FO(3) = CXZERO
-                     FO(4) = DCMPLX( OMEGAM, RXZERO )
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS - P(3) )
-                     FO(1) = OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FO(2) = OMEGAP * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( P(1), -P(2) )
-                     FO(3) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FO(4) = OMEGAM * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( P(1), -P(2) )
-                  END IF
-               END IF
-            ELSE IF ( NHEL .EQ. -1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FO(1) = CXZERO
-                     FO(2) = DCMPLX( OMEGAM, RXZERO )
-                     FO(3) = CXZERO
-                     FO(4) = DCMPLX( OMEGAP, RXZERO )
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS + P(3) )
-                     FO(1) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( -P(1), -P(2) )
-                     FO(2) = OMEGAM * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FO(3) = OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( -P(1), -P(2) )
-                     FO(4) = OMEGAP * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FO(1) = DCMPLX( -OMEGAM, RXZERO )
-                     FO(2) = CXZERO
-                     FO(3) = DCMPLX( -OMEGAP, RXZERO )
-                     FO(4) = CXZERO
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS - P(3) )
-                     FO(1) = OMEGAM * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( -P(1), -P(2) )
-                     FO(2) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FO(3) = OMEGAP * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( -P(1), -P(2) )
-                     FO(4) = OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                  END IF
-               END IF
-            ELSE
-               STOP 'OXXXXX:  FERMION HELICITY MUST BE +1,-1'
-            END IF
-         ELSE IF ( NSF .EQ. -1 ) THEN
-            IF ( NHEL .EQ. 1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FO(1) = CXZERO
-                     FO(2) = DCMPLX( OMEGAM, RXZERO )
-                     FO(3) = CXZERO
-                     FO(4) = DCMPLX( -OMEGAP, RXZERO )
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS + P(3) )
-                     FO(1) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( -P(1), -P(2) )
-                     FO(2) = OMEGAM * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FO(3) = -OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( -P(1), -P(2) )
-                     FO(4) = -OMEGAP * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FO(1) = DCMPLX( -OMEGAM, RXZERO )
-                     FO(2) = CXZERO
-                     FO(3) = DCMPLX( OMEGAP, RXZERO )
-                     FO(4) = CXZERO
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS - P(3) )
-                     FO(1) = OMEGAM * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( -P(1), -P(2) )
-                     FO(2) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FO(3) = -OMEGAP * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( -P(1), -P(2) )
-                     FO(4) = -OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                  END IF
-               END IF
-            ELSE IF ( NHEL .EQ. -1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FO(1) = DCMPLX( -OMEGAP, RXZERO )
-                     FO(2) = CXZERO
-                     FO(3) = DCMPLX( OMEGAM, RXZERO )
-                     FO(4) = CXZERO
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS + P(3) )
-                     FO(1) = -OMEGAP * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FO(2) = -OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( P(1), -P(2) )
-                     FO(3) = OMEGAM * RS2PA
-     &                     * DCMPLX( SPAZ, RXZERO )
-                     FO(4) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( P(1), -P(2) )
-                  END IF
-               ELSE
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FO(1) = CXZERO
-                     FO(2) = DCMPLX( -OMEGAP, RXZERO )
-                     FO(3) = CXZERO
-                     FO(4) = DCMPLX( OMEGAM, RXZERO )
-                  ELSE
-                     RS2PA = RXONE / SQRT( RXTWO * PABS )
-                     SPAZ = SQRT( PABS - P(3) )
-                     FO(1) = -OMEGAP * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FO(2) = -OMEGAP * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( P(1), -P(2) )
-                     FO(3) = OMEGAM * RS2PA / SPAZ
-     &                     * DCMPLX( PLAT, RXZERO )
-                     FO(4) = OMEGAM * RS2PA * SPAZ / PLAT
-     &                     * DCMPLX( P(1), -P(2) )
-                  END IF
-               END IF
-            ELSE
-               STOP 'OXXXXX:  FERMION HELICITY MUST BE +1,-1'
-            END IF
-         ELSE
-            STOP 'OXXXXX:  FERMION TYPE MUST BE +1,-1'
-         END IF
-C
-C do massless case
-C
-      ELSE
-         IF ( NSF .EQ. 1 ) THEN
-            IF ( NHEL .EQ. 1 ) THEN
-               IF ( P(3) .GE. RXZERO ) THEN
-                  IF ( PLAT .EQ. RXZERO ) THEN
-                     FO(1) = DCMPLX( OMEGAP, RXZERO )
-                     FO(2) = CXZERO
-                     FO(3) = CXZERO
-                     FO(4) = CXZERO
-                  ELSE
-                     SPAZ = SQRT( PABS + P(3) )
-                     FO(1) = DCMPLX( SPAZ, RXZERO )
-                     FO(2) = RXONE / SPAZ
-     &                     * DCMPLX( P(1), -P(2) )
-