#!/bin/bash
# finds out architecture name, defines system files, and construct
# the relevant directory for run, if not present. If the architecture
# cannot be figured out, defines ./Run as working directory, and use 
# sun as system file -- it should work in most of the cases
function architecture {
PATH=".:""$PATH"
thisdir=`pwd`
homedir=$HOME
case `uname` in
Linux*) ARCHNAME=LinuxPY8 ; DBGFLAG='DEBUG=-ggdb' ; EXTRA='EXTRAOBJ=linux.o' ;
        EXMAKE=make ; CPOPT=pu ;;
OSF*)  ARCHNAME=AlphaPY8 ; DBGFLAG='DEBUG=-g' ; EXTRA='EXTRAOBJ=alpha.o' ;
       EXMAKE=gmake ; CPOPT=p ;;
Sun*) ARCHNAME=SunPY8 ; DBGFLAG='DEBUG=-g' ; EXTRA='EXTRAOBJ=sun.o' ;
      EXMAKE=gmake ; CPOPT=p ;;
Darwin*) ARCHNAME=DarwinPY8 ; DBGFLAG='DEBUG=-g' ; EXTRA='EXTRAOBJ=sun.o' ;
         EXMAKE=make ; CPOPT=p ;;
*) ARCHNAME=RunPY8 ; DBGFLAG='DEBUG=-g' ; EXTRA='EXTRAOBJ=sun.o' ;;
esac
#
}

# $? is the value of last executed command. A call to this function
# after a failure will cause the program to quit the script
function teststatus {
rc=$?
if [ 0 = $rc ]
then
:
else
echo $* did not succeed, exit status=$rc 
exit $rc
fi
}

# returns a string which identifies the target
# for the Makefile
function getSTRmake {
stringmake=MadFKS
}

# utility function for dothelinks
function stripextension {
echo $1 | sed "s/\..*\$//"
}

# utility function for dothelinks
function capitalize {
echo $1 | sed "y/abcdefghijklmnopqrstuvwxyz/ABCDEFGHIJKLMNOPQRSTUVWXYZ/"
}

# creates logical links for the PDF grid files. By P. Nason
function dothelinks {
if [ -d $PDFPATH ]
then
for i in ${PDFPATH}/*.dat ${PDFPATH}/*.tbl
do
if [ -f $i ]
then
name=`basename $i`
name=`stripextension $name`
case $name in
mrst200*) ;;
*mrs*|grpol|grvdm|lac1|pion[1-3]) name=`capitalize $name`;;
esac
if [ ! -L $thisdir/$name ] || [ ! $thisdir/$name -ef $i ]
then
ln -sf $i $thisdir/$name
fi
fi
done
for i in ${PDFPATH}/a02.*
do
if [ -f $i ]
then
name=`basename $i`
if [ ! -L $thisdir/$name ] || [ ! $thisdir/$name -ef $i ]
then
ln -sf $i $thisdir/$name
fi
fi
done
fi
}

# creates logical links for LHAPDF, and replaced PDF group name (unused
# by LHAPDF) with a LHAPDF-specific string
function linklhapdf {
case $LHAOFL in
FREEZE|freeze) PDFGROUP=LHAPDF ;;
EXTRAPOLATE|extrapolate) PDFGROUP=LHAEXT ;;
*) echo "no such option; failure in linklhapdf" ; exit 1 ;;
esac
}    

# creates the executable for NLO and MC; the entry of the function
# is either NLO or MC, depending on which executable one wants to obtain;
# PDFLIB or private PDF library is linked, depending on the value of
# the parameter UsedPdfLib, which is set by the function whichpdflib.
# Extra libraries are listed in the shell variable EXTRALIBS; do NOT
# remove the double quotes in the call to nameextralib, otherwise only
# the first library of the list is taken
function compile {
architecture

EVTDIR=`pwd`"/../Events"

if [ "$PDFCODE" -ne 0 ]
then
  PYTPDF='EXTPDF'
  PDFLIBRARY='LHAPDF'
  UsedPdfLib='LHAPDF'
  export LIBLHAPDF=$LHAPDFPATH/lib/libLHAPDF.a
  LHALINK=DYNAMIC
  if [ "$PDFCODE" -gt 0 ]; then LHAOFL=FREEZE; fi
  if [ "$PDFCODE" -gt 1 ]; then PDFSET=$PDFCODE; fi
  if [ "$PDFCODE" -lt 0 ]; then LHAOFL=EXTRAPOLATE; fi
  if [ "$PDFCODE" -lt -1 ]; then PDFSET=$((PDFCODE*-1)); fi
  linklhapdf
# convert PDF number to name reading PDFsets.index
  ISLHA5=0
  ISLHA6=0
  LHAVERSION=`$LHAPDFPATH/bin/lhapdf-config --version`
  if [ "${LHAVERSION:0:1}" == "5" ]
  then
      LHAINDEX=$LHAPDFPATH/share/lhapdf/PDFsets.index
      ISLHA5=1
  elif [ "${LHAVERSION:0:1}" == "6" ]
  then
      LHAINDEX=$LHAPDFPATH/share/LHAPDF/pdfsets.index
      ISLHA6=1
  else
      echo 'Unable to understand LHAPDF version!'
      exit 1
  fi
  while read pdfline
  do
    pdfinfo=($(echo ${pdfline}))
    if [ "$PDFSET" == "${pdfinfo[0]}" ]
    then
	if [ "$ISLHA5" -eq 1 ]; then PDFSETN="${pdfinfo[4]}"; fi
	if [ "$ISLHA6" -eq 1 ]; then PDFSETN="${pdfinfo[1]}"; fi
    fi
  done < $LHAINDEX
elif [ "$PDFCODE" -eq 0 ]
then
  PYTPDF='DEFAULT'
  PDFLIBRARY='THISLIB'
  UsedPdfLib='THISLIB'
  export LIBLHAPDF=-llhapdfdummy
# the following is dummy
  LHALINK=DYNAMIC
  LHAOFL=FREEZE
else
  echo 'Invalid PDFCODE' $PDFCODE
  exit 1
fi

whichpdflib $PDFLIBRARY
whichlhapdf $LHALINK
# stdhep is not used: do not look for libraries
IFS=' ' read -a exlibs <<< "${EXTRALIBS}"
for i in ${!exlibs[@]}
do
 if [[ "${exlibs[i]}" == "stdhep" || "${exlibs[i]}" == "Fmcfio" ]]
 then
   unset exlibs[$i]
 fi
done
EXTRALIBS="$(echo ${exlibs[@]})"
nameextralib "$EXTRALIBS"
pathextralib "$EXTRAPATHS"
pathextralibdyn "$EXTRAPATHS"
pathextraincl "$INCLUDEPATHS"
EXEC_NAME=
EXEC_NAMEF=
EXEC_NAMEFF=
PY8UTIMAKE="PY8UTI="$PY8UTI
PY8MAKE="PY8VER="$PY8VER
LIBSMAKE="EXTRALIBS="$ExtraLibs
LIBSMAKEP="EXTRAPATHS="$ExtraPaths
INCLMAKE="INCLOPTION="$ExtraIncl
INCDIRMK="INCDIR="`pwd`"/include"
SRCDIRMK="SRCDIR="`pwd`"/srcPythia8"
COMSRCMK="COMSRC="`pwd`"/srcCommon"
EVTDIR=`pwd`"/../Events"
if [ "$LD_LIBRARY_PATH" = "" ]
then
  LD_LIBRARY_PATH="$ExtraPathsDyn"
else
  if [ "$ExtraPathsDyn" != "" ]
  then
    LD_LIBRARY_PATH="$ExtraPathsDyn"":""$LD_LIBRARY_PATH"
  fi
fi
LIBSLHA="LHALIB="
case $1 in
NLO) echo "NLO must be computed with MadFKS"; exit 1 ;;
MC) if [ $UsedPdfLib = "PDFLIB" ]
     then
       EXEC_NAMEF="PY8_EXE_DEFAULT"
     elif [ $UsedPdfLib = "THISLIB" ]
     then
       EXEC_NAMEF="PY8_EXE_DEFAULT"
     elif [ $UsedPdfLib = "LHAPDF" ]
     then
       EXEC_NAMEF="PY8_EXE_LHAPDF"
     else
       echo "fatal error in compile"; exit 1
     fi ;;
*) echo "do not know what to do in function compile" ; exit 1 ;;
esac
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HEPMCPATH/lib
if [ $1 = "NLO" ]
then
  getSTRmake
  EXEC_NAME=$stringmake$EXEC_NAMEF
elif [ $1 = "MC" ]
then
  EXEC_NAME="$EXEC_NAMEF"
fi
EXEC_NAMEFF="$EXEPREFIX""$EXEC_NAMEF"

rm -f $thisdir/$EXEC_NAMEFF
# new compilation step for PY8 mode
export PYTHIA8LOCATION=$PY8PATH

# check the PYTHIA8 version by the position of certain files

if [ ! -d $PYTHIA8LOCATION/xmldoc ]
then
 PY8VER="8.2"
 echo ' USING PYTHIA 8, version >= 8.200'
 chmod +x $PYTHIA8LOCATION/bin/pythia8-config
 read HEPMCINCLIB <<< $($PYTHIA8LOCATION/bin/pythia8-config --hepmc2)
else
 PY8VER="8.1"
 echo ' USING PYTHIA 8, version < 8.200'
fi

if [ "$PY8UTI" != "" ]
then
 if [ $PY8VER == "8.1" ]; then cp -rf $thisdir/srcPythia8/Pythia8.cc $thisdir/Pythia8.cc; fi
 if [ $PY8VER == "8.2" ]; then cp -rf $thisdir/srcPythia8/Pythia82.cc $thisdir/Pythia82.cc; fi
 cp -rf $thisdir/srcPythia8/Makefile $thisdir/Makefile
else
 if [ $PY8VER == "8.1" ]; then cp -rf $thisdir/srcPythia8/Pythia8_hep.cc $thisdir/Pythia8.cc; fi
 if [ $PY8VER == "8.2" ]; then cp -rf $thisdir/srcPythia8/Pythia82_hep.cc $thisdir/Pythia82.cc; fi
 cp -rf $thisdir/srcPythia8/Makefile_hep $thisdir/Makefile
fi
if [ -f $thisdir/Pythia8.exe ]
then
  rm -rf $thisdir/Pythia8.exe
fi

if [ $PY8VER == "8.1" ]; then cp -rf $PYTHIA8LOCATION/xmldoc $thisdir; fi
if [ $PY8VER == "8.2" ]; then cp -rf $PYTHIA8LOCATION/share/Pythia8/xmldoc $thisdir; fi
PY8UTIMAKE="PY8UTI="$PY8UTI
FORFILEMAKE="FORFILE="$PY8UTI
HEPMCMAKE="HEPMCINCLIB="$HEPMCINCLIB
LGFORTRANMAKE="LIBGFORTRAN=-lgfortran"

if [ $PY8VER == "8.1" ]; then make Pythia8 "$PY8UTIMAKE" "$FORFILEMAKE" "$LGFORTRANMAKE" "$LIBSMAKE" "$LIBSMAKEP" "$INCLMAKE"; fi
if [ $PY8VER == "8.2" ]; then make Pythia82 "$PY8UTIMAKE" "$FORFILEMAKE" "$LGFORTRANMAKE" "$LIBSMAKE" "$LIBSMAKEP" "$INCLMAKE" "$HEPMCMAKE"; fi

cd $thisdir
if [ ! -f $thisdir/Pythia8.exe ]
then
  echo "Pythia8 compilation did not succeed, exiting"
  exit -1
fi


teststatus Compilation
}

# foolproof utility function, to be called by the user
function compileMC {
compile MC
}

# wrapper for the functions which run and compile MC codes
function runMC {
compile MC

#put MCMODE to upper case
MCMODE=`echo $MCMODE | tr '[a-z]' '[A-Z]' `
# Check correspondence between MCMODE and what written in the event file
if [ $MCMODE != "PYTHIA8" ]
then
  echo "Error! PYTHIA8 Script called with MCMODE not equal to PYTHIA8"
  exit 1
fi

UMASSPY=$UMASS
DMASSPY=$DMASS
SMASSPY=$SMASS
CMASSPY=$CMASS
if (( $(bc <<< "$B_MASS <  0.0") )); then BMASSPY=$BMASS; fi
if (( $(bc <<< "$B_MASS >= 0.0") )); then BMASSPY=$B_MASS; fi
EMASSPY=$EMASS
MUMASSPY=$MUMASS
TAUMASSPY=$TAUMASS
GMASSPY=$GMASS

if [ $BEAM1 == 1 ]; then iA=2212; elif [ $BEAM1 == -1 ]; then iA=-2212;
elif [ $BEAM1 == 0 ]; then iA=-11; elif [ $BEAM1 == 2 ]; then iA=2112;
elif [ $BEAM1 == -2 ]; then iA=-2112; else echo "Unrecognized incoming beam 1";
exit 1; fi
if [ $BEAM2 == 1 ]; then iB=2212; elif [ $BEAM2 == -1 ]; then iB=-2212;
elif [ $BEAM2 == 0 ]; then iB=11; elif [ $BEAM2 == 2 ]; then iB=2112;
elif [ $BEAM2 == -2 ]; then iB=-2112; else echo "Unrecognized incoming beam 2";
exit 1; fi
# Check the previous!

whichpytpdf $PYTPDF
if [ $pdftype = 1 ]
then
  whichpdflib $PDFLIBRARY
  if [ $UsedPdfLib = "THISLIB" ]
  then
    dothelinks
  elif [ $UsedPdfLib = "LHAPDF" ]
  then
    linklhapdf
  fi
fi
ERR_NUM_PY8=$(echo "$ERR_FR_PY8 * $NEVENTS" | bc)
ERR_NUM_PY8=${ERR_NUM_PY8/.*}
let "ERR_NUM_PY8 += 1"

echo " "
echo "*****   Now ready for showering" $NEVENTS "events with $MCMODE "
echo " "
runMCMadFKS
teststatus runMC
}

# compiles and runs the MC
function runMCMadFKS {
ifile="MCATNLO_"$MCMODE"_input"
exefile="MCATNLO_"$MCMODE"_EXE"
if [ -f ./$ifile ]
then
\rm ./$ifile 
fi
#mv $EXEC_NAMEFF $exefile
# create input file
cat <<EOF > $ifile

! 1) Settings used in the main program.
Main:numberOfEvents = $NEVENTS_TOT    ! Number of events in the LHE file
Main:spareMode1 = $NEVENTS            ! Number of events to be showered
Main:spareWord1 = $EVENT_NORM         ! Event weights are normalized to sum 
                                      ! or average to the cross section
Main:timesAllowErrors = $ERR_NUM_PY8  ! Number of allowed errors

! 2) Settings related to output in init(), next(), amd stat().
Init:showChangedParticleData = off    ! Shows all non-default particle settings
Next:numberShowEvent = $MAXPR_PY8     ! print event record n times
Stat:showProcessLevel = on            ! Process statistics
Stat:showErrors = on                  ! Error statistics
Check:epTolErr = 0.001                ! Momentum-conservation tolerance 

! 3) Beam-parameter settings.
Beams:frameType = 4                ! LHE initialization
Beams:LHEF = events.lhe            ! Input LHE file

! 4) Switch on/off the key event-generation steps.
EOF
if [ $UsedPdfLib = "LHAPDF" ]
then
 if [ $PY8VER == "8.1" ]
 then
  cat <<EOF >> $ifile
PDF:useLHAPDF = on                 ! Use of LHAPDF
PDF:LHAPDFset = $PDFSETN           ! PDF set
EOF
  if [ $PDFGROUP = "LHAPDF" ]
  then
  cat <<EOF >> $ifile
PDF:extrapolateLHAPDF = off        ! extrapolate PDF set outside the boundaries
EOF
  elif [ $PDFGROUP = "LHAEXT" ]
  then
  cat <<EOF >> $ifile
PDF:extrapolateLHAPDF = on         ! extrapolate PDF set outside the boundaries
EOF
  else
  echo "Unknown PDFGROUP " $PDFGROUP
  exit 1
  fi
 else
  if [ "$ISLHA5" -eq 1 ]
  then
  cat <<EOF >> $ifile
PDF:pSet = LHAPDF5:$PDFSETN    ! PDF set used in showers and multiparton interactions
EOF
   if [ $PDFGROUP = "LHAPDF" ]
   then
   cat <<EOF >> $ifile
PDF:extrapolate = off        ! extrapolate PDF set outside the boundaries
EOF
   elif [ $PDFGROUP = "LHAEXT" ]
   then
   cat <<EOF >> $ifile
PDF:extrapolate = on         ! extrapolate PDF set outside the boundaries
EOF
   else
   echo "Unknown PDFGROUP " $PDFGROUP
   exit 1
   fi
  elif [ "$ISLHA6" -eq 1 ]
  then
  cat <<EOF >> $ifile
PDF:pSet = LHAPDF6:$PDFSETN    ! PDF set used in showers and multiparton interactions
EOF
  fi
 fi
else
 if [ $PY8VER == "8.1" ]
 then
 cat <<EOF >> $ifile
PDF:useLHAPDF = off                ! Use of LHAPDF
EOF
 fi
 cat <<EOF >> $ifile
PDF:pSet = 13                      ! PDF set used in showers and multiparton interactions (13 = NNPDF2.3 QCD+QED)
EOF
fi

if ((`bc <<< "$LAMBDAPYTH >= 0.0"`))
then
cat <<EOF >> $ifile
#coupSM:Lambda5 = $LAMBDAPYTH      ! five-flavour LambdaQCD
EOF
fi

if [ $UE_PY8 = ".FALSE." ]
then
cat <<EOF >> $ifile
PartonLevel:MPI = off              ! Multiparton interactions
                                   ! ** USE on FOR REALISTIC SHOWER SIMULATION ** !!!
EOF
else
cat <<EOF >> $ifile
PartonLevel:MPI = on               ! Multiparton interactions
EOF
fi

if [ $HADRONIZE_PY8 = ".FALSE." ]
then
cat <<EOF >> $ifile
HadronLevel:Hadronize = off         ! Hadronisation
                                   ! ** USE on FOR REALISTIC SHOWER SIMULATION ** !!!
EOF
else
cat <<EOF >> $ifile
HadronLevel:Hadronize = on         ! Hadronisation
EOF
fi

cat <<EOF >> $ifile

! 5) Final-state shower.
TimeShower:pTmaxMatch = 1          ! Use scalup (re-check)
TimeShower:pTmaxFudge = 1.         ! Factor changing the max scale
TimeShower:alphaSvalue = 0.118     ! Alpha_s(MZ) in final-state shower
TimeShower:alphaSorder = 1         ! Alpha_s running order in final-state shower
TimeShower:alphaEMorder = 0        ! Alpha_EM running order in final-state shower
TimeShower:dampenBeamRecoil = off  ! Dampens the effect of beam recoil
TimeShower:globalRecoil = on       ! All final-state particles recoil against the branching
TimeShower:nMaxGlobalRecoil = 1    ! Number of splittings with TimeShower:globalRecoil = on
TimeShower:globalRecoilMode = 2    ! Global recoil only for S events whose first emission is FSR
TimeShower:nMaxGlobalBranch = 1    ! Number of FSR splittings proposed with global recoil
TimeShower:nPartonsInBorn = -1      ! Number of Born QCD final-state partons (to treat H and S differently)
TimeShower:limitPTmaxGlobal = on   ! Limits pT < min(SCALUP,mDipole/2).

TimeShower:QEDshowerByQ = off      ! Prevent quarks from radiating photons
                                   ! ** USE on FOR REALISTIC SHOWER SIMULATION ** !!!
TimeShower:QEDshowerByL = off      ! Prevent leptons from radiating photons
                                   ! ** USE on FOR REALISTIC SHOWER SIMULATION ** !!!
TimeShower:QEDshowerByGamma = off  ! Prevent photons from branching
                                   ! ** USE on FOR REALISTIC SHOWER SIMULATION ** !!!

TimeShower:MEcorrections = off     ! No Matrix-element corrections
TimeShower:alphaSuseCMW = false    ! Use the CMW prescription in FSR
TimeShower:weightGluonToQuark = 1  ! Use normal Altarelli-Parisi kernels for g -> q qbar

! 6) Initial-state shower.
SpaceShower:pTmaxMatch = 1         ! Use scalup (re-check)
SpaceShower:pTmaxFudge = 1.        ! Factor changing the max scale
SpaceShower:alphaSvalue = 0.118    ! Alpha_s(MZ) in final-state shower
SpaceShower:alphaSorder = 1        ! Alpha_s running order in final-state shower
SpaceShower:alphaEMorder = 0       ! Alpha_EM running order in final-state shower

SpaceShower:QEDshowerByQ = off     ! Prevent quarks from radiating photons
                                   ! ** USE on FOR REALISTIC SHOWER SIMULATION ** !!!
SpaceShower:QEDshowerByL = off     ! Prevent leptons from radiating photons
                                   ! ** USE on FOR REALISTIC SHOWER SIMULATION ** !!!

SpaceShower:MEcorrections = off    ! No Matrix-element corrections
SpaceShower:rapidityorder = off    ! Do not order branchings in rapidity
SpaceShower:alphaSuseCMW = false   ! Use the CMW prescription in ISR

! 7) Non-perturbative stuff.
BeamRemnants:primordialKT = off    ! No primordial kT
                                   ! ** USE on FOR REALISTIC SHOWER SIMULATION ** !!!

! 8) Particle features.
1:m0 = $DMASSPY                    ! down mass
2:m0 = $UMASSPY                    ! up mass
3:m0 = $SMASSPY                    ! strange mass
4:m0 = $CMASSPY                    ! charm mass
5:m0 = $BMASSPY                    ! bottom mass
6:m0 = $TMASS                      ! top mass
11:m0 = $EMASSPY                   ! electron mass
13:m0 = $MUMASSPY                  ! muon mass
15:m0 = $TAUMASSPY                 ! tauon mass
23:m0 = $ZMASS                     ! Z mass
24:m0 = $WMASS                     ! W mass
25:m0 = $HGGMASS                   ! Higgs mass
6:mWidth = $TWIDTH                 ! top width
23:mWidth = $ZWIDTH                ! Z width
24:mWidth = $WWIDTH                ! W width
25:mWidth = $HGGWIDTH              ! Higgs width
EOF

if [ $ICKKW = 0 ]
then
cat <<EOF >> $ifile

! 9) Multijet merging.
JetMatching:doFxFx = off          ! switch off FxFx


EOF
elif [ $ICKKW = 3 ]
then
cat <<EOF >> $ifile

! 9) Multijet merging.
JetMatching:doFxFx = on           ! switch on FxFx
JetMatching:merge = on            ! switch on MLM-merging machinery
JetMatching:qCut = $QCUT          ! merging scale
JetMatching:qCutME = $PTJCUT      ! minimum-jet(pT) cut
JetMatching:coneRadius = 1.0      ! jet radius definition
JetMatching:etaJetMax = 1000.0    ! max jet rapidity definition
JetMatching:nJetMax = $NJMAX      ! max multiplicity
JetMatching:scheme = 1            ! MadGraph MLM-type merging
JetMatching:setMad = off          ! Don't read from LHE in MG5LO format

EOF
fi
cat <<EOF >> $ifile
! 10) Decays.
EOF

if [ $PI_STABLE_PY8 = ".TRUE." ]
then
cat <<EOF >> $ifile
111:mayDecay = false               ! stable pi0
EOF
fi
if [ $B_STABLE_PY8 = ".TRUE." ]
then
cat <<EOF >> $ifile
511:maydecay = false               ! stable B hadrons 
521:maydecay = false               ! stable B hadrons
531:maydecay = false               ! stable B hadrons
541:maydecay = false               ! stable B hadrons
551:maydecay = false               ! stable B hadrons
553:maydecay = false               ! stable B hadrons
5112:maydecay = false              ! stable B hadrons
5122:maydecay = false              ! stable B hadrons
5132:maydecay = false              ! stable B hadrons
5222:maydecay = false              ! stable B hadrons
5232:maydecay = false              ! stable B hadrons
5332:maydecay = false              ! stable B hadrons
EOF
fi
if [[ $WP_STABLE_PY8 = ".TRUE." || $WM_STABLE_PY8 = ".TRUE." ]]
then
cat <<EOF >> $ifile
24:maydecay = false                ! stable W boson
EOF
fi
if [ $Z_STABLE_PY8 = ".TRUE." ]
then
cat <<EOF >> $ifile
23:maydecay = false                ! stable Z boson
EOF
fi
if [ $H_STABLE_PY8 = ".TRUE." ]
then
cat <<EOF >> $ifile
25:maydecay = false                ! stable Higgs boson
EOF
fi
if [ $TAUP_STABLE_PY8 = ".TRUE." ]
then
cat <<EOF >> $ifile
-15:maydecay = false               ! stable tau+
EOF
fi
if [ $TAUM_STABLE_PY8 = ".TRUE." ]
then
cat <<EOF >> $ifile
15:maydecay = false                ! stable tau-
EOF
fi
if [ $MUP_STABLE_PY8 = ".TRUE." ]
then
cat <<EOF >> $ifile
-13:maydecay = false               ! stable mu+
EOF
fi
if [ $MUM_STABLE_PY8 = ".TRUE." ]
then
cat <<EOF >> $ifile
13:maydecay = false                ! stable mu-
EOF
fi
if [ $RNDEVSEED_PY8 != 0 ]
then
cat <<EOF >> $ifile
Random:setSeed = on                ! random seed
Random:seed = $RNDEVSEED_PY8
EOF
fi


arr=()
numDM=0
iDM=0
for i in {1..99}
do
 eval arr_elem='$'DM_$i
 if [ "$arr_elem" != "" ]
 then
  arr=("${arr[@]}" "$arr_elem")
 fi
done

for DM in "${arr[@]}"
do
cat <<EOF >> ./$ifile
$DM  ! decay mode
EOF
done


cp -f $thisdir/$ifile Pythia8.cmd

}

# this function set the parameter pdftype according to the value 
# of PYTPDF (the entry of this function) given in input
function whichpytpdf {
case $1 in
DEFAULT|default) pdftype=0 ;;
EXTPDF|extpdf) pdftype=1 ;;
*) echo "error in whichpytpdf: no such option" ; exit 1 ;;
esac
}

# checks that the value given to PDFLIBRAY in input is meaningful
function whichpdflib {
case $1 in
THISLIB|thislib) UsedPdfLib=THISLIB ;;
PDFLIB|pdflib) UsedPdfLib=PDFLIB ;;
LHAPDF|lhapdf) UsedPdfLib=LHAPDF ;;
*) echo "no such library for PDFS; failure in whichpdflib" ; exit 1 ;;
esac
}

# checks that the value given to LHALINK in input is meaningful
function whichlhapdf {
case $1 in
STATIC|static) UsedLhaPdf=lhasta ;;
DYNAMIC|dynamic) UsedLhaPdf=lhadyn ;;
*) echo "no such option for LHAPDF; failure in whichlhapdf" ; exit 1 ;;
esac
}

# prepends -l to library names
function nameextralib {
ExtraLibs="$1"
ilen1=${#ExtraLibs}
if [ $ilen1 -ne 0 ]
then
  ExtraLibs=`echo $ExtraLibs | sed "s/ / -l/g"`
  ExtraLibs="-l""$ExtraLibs"
fi
}

# prepends -L to library paths
function pathextralib {
ExtraPaths="$1"
ilen2=${#ExtraPaths}
if [ $ilen2 -ne 0 ]
then
  ExtraPaths=`echo $ExtraPaths | sed "s/ / -L/g"`
  ExtraPaths="-L""$ExtraPaths"
fi
}

# separate library names with columns
function pathextralibdyn {
ExtraPathsDyn="$1"
ilen3=${#ExtraPathsDyn}
if [ $ilen3 -ne 0 ]
then
  ExtraPathsDyn=`echo $ExtraPathsDyn | sed "s/ /:/g"`
fi
}

# prepends -I to include paths
function pathextraincl {
ExtraIncl="$1"
ilen4=${#ExtraIncl}
if [ $ilen4 -ne 0 ]
then
  ExtraIncl=`echo $ExtraIncl | sed "s/ / -I/g"`
  ExtraIncl="-I""$ExtraIncl"
fi
}


# Collection of Pythia8 settings that used to appear in the
# card, and that are supefluous because overridden or default

#Beams:idA = $iA                    ! Beam identities
#Beams:idB = $iB                    ! Beam identities
#ProcessLevel:all = on              ! Generation
#ProcessLevel:resonanceDecays = on  ! Resonance decays
#PartonLevel:all = on               ! Paron level: if off, stops after hard
#                                   ! process generation
#PartonLevel:ISR = on               ! Initial state shower
#PartonLevel:FSR = on               ! Final state shower
#PartonLevel:FSRinProcess = on      ! Final state shower in association with
#                                   ! the hard process
#PartonLevel:FSRinResonances = on   ! Final state shower in resonance decays
#HadronLevel:all = on               ! Hadron level: if off, stops before
#                                   ! hadronization
#HadronLevel:Decay = on             ! Hadron decays
#PhaseSpace:mHatMin = 4.            ! Min invariant mass
#PhaseSpace:mHatMax = -1.           ! Max invariant mass
#PhaseSpace:pTHatMin = 0.           ! Min pT in 2->2
#PhaseSpace:pTHatMax = -1.          ! Max pT in 2->2
#PhaseSpace:pTHatMinDiverge = 1.    ! If massless final state, to avoid
#                                   ! divergences
#PhaseSpace:useBreitWigners = on    ! Masses according to Breit-Wigner
#PhaseSpace:pTHat3Min = 0.          ! Min pT for the hardest parton in 2->3
#PhaseSpace:pTHat3Max = -1.         ! Max pT for the hardest parton in 2->3
#PhaseSpace:pTHat5Min = 0.          ! Min pT for the softest parton in 2->3
#PhaseSpace:pTHat5Max = -1.         ! Max pT for the softest parton in 2->3
#PhaseSpace:RsepMin = 0.            ! Min R separation in 2->3
#TimeShower:interleave = on         ! If on, FSR interleaved with ISR
#TimeShower:allowBeamRecoil = on    ! If off, no energy transfer to from
#                                   ! ISR to FSR
#TimeShower:QCDshower = on          ! QCD final-state shower
#TimeShower:nGluonToQuark = 5       ! Number if flavors allowed in g->qqbar
#TimeShower:phiPolAsym = on         ! Azimuthal asymmetry induced by gluon
#                                   ! polarization
#SpaceShower:QCDshower = on         ! QCD final-state shower
#SpaceShower:phiPolAsym = on        ! Azimuthal asymmetry induced by gluon
#                                   ! polarization
#SpaceShower:nQuarkIn = 5           ! Number of flavors in g->qqbar and also in
#                                   ! incoming beams
#Init:showChangedSettings = on      ! Shows all non-default settings
#Next:numberCount = 100             ! print message every n events
#Next:numberShowInfo = 2            ! print event information n times
#Next:numberShowProcess = 1         ! print process record n times
#TimeShower:MEafterFirst = off      ! No Matrix-element corrections after first emission
#                                   ! it is already off if MEcorrections = off
#SpaceShower:MEafterFirst = off     ! No Matrix-element corrections after first emission
#                                   ! it is already off if MEcorrections = off