#!/bin/bash

################################################################################
# General settings for digitization/reconstruction                             #
################################################################################
iNLogicalCores=1           # number of threads used to compute the charge memory coefficients
                           # 1 core reserved for the main thread
iTimeLimit=0-05:00:00      # expected digitization job completion time (TIMEOUT if exceeded)
cPartition=long            # Kronos partitions (long, main, debug)

cSetup=sps_feb15_sim       # as defined in geometry/setup/setup_${cSetup}.C
cGeometry=v15i             # needed for documentation purposes only
                           # should match the string 'tofGeoTag' in the setup file
iNParallelSpills=10        # number of parallel digitization jobs
iNSerialSpills=1           # number of consecutive spills in a single digitization job
bAverageSerialSpills=1     # 1: clear MRPC memory on every new spill; time in spill
                           # 0: never clear MRPC memory; time in run
dSpillLength=10.           # [s]
dSpillBreakLength=20.      # [s] must be significantly larger than the mean time between events
bSingleTimeSlice=0         # WARNING: single ROOT branch entries CANNOT contain more than 1 GiB of data!!!
                           #          enable single TS mode only if you are sure that your 'TClonesArray'
                           #          buffers holding (unsplittable) 'CbmMatch' objects will not exceed
                           #          the 1 GiB limit.
cRunMode=hic               # scf: single-counter flooding
                           #      illuminate the specified counter with a box source
                           # haf: homogeneous area flooding
                           #      illuminate the (entire) setup with a box source
                           # hic: heavy-ion collision
                           #      illuminate the (entire) setup with HI collisions
bRunWithMemory=0           # 1: memorize all ToF points generated in a counter and
                           #    calculate during digitization for each point processed
                           #    the spatiotemporal impact of all preceding points
                           #    on the local detector response function
                           # 0: do not memorize any ToF point, i.e. the detector
                           #    response will not degrade as a function of irradiation
                           #    time

################################################################################
# Default common run parameters (to be overwritten by each run mode)           #
################################################################################
bRunReco=1                 # 1: run the reconstruction macro
bRunQA=1                   # 1: add 'CbmTofRecoMCQa' to the list of reconstruction tasks
                           #    and process QA histograms in 'combine_histograms.sh'
bAlternativeBranchNames=0  # 1: prepend an "A" to all data branch names to avoid
                           #    equally named input and output branches during
                           #    event building at the beginning of the analysis
                           #    chain
bBuildCbmEvents=0          # 1: add 'CbmTofBuildMCEvents' to the list of reconstruction tasks
bBuildHits=1               # 1: add 'CbmTofClusterMC' and 'CbmTofMatchReco' to the list of
                           #    reconstruction tasks
bHitCompatibility=0        # 1: hit building compatibility with 'CbmTofTestBeamClusterizer'
bSpotResponse=0            # 1: spot illumination
                           # 0: flood illumination (always in "hic" mode)
bNoTransportPhysics=0      # 1: disable all physics processes during MC transport
                           # 0: perform a regular MC simulation
bTwoSidedReadout=0         # 1: enforces a two-sided readout for debugging purposes even
                           #    of counters that are explicitly declared one-sided in the
                           #    corresponding *.digitb.par file
bCalibrateDigis=0          # 1: write ideally calibrated digi objects to output (time and ToT)
bStorePointsInTS=0         # 1: fill MC points into time slice output
                           #    (requires a double precision 'FairMCPoint::fTime' !!!)

################################################################################
# SCF-specific settings                                                        #
################################################################################
if [ "${cRunMode}" == "scf" ]; then
  bSingleCounter=1
  dParticleFlux=30000.     # [Hz/cm^2]
#  iModuleType=3
#  iModuleIndex=0
#  iCounterIndex=0
#  dGlassDimX=32.           # [cm]
#  dGlassDimY=26.9          # [cm]
#  iModuleType=4
#  iModuleIndex=0
#  iCounterIndex=0
#  dGlassDimX=15.           # [cm]
#  dGlassDimY=4.            # [cm]
#  iModuleType=5
#  iModuleIndex=2
#  iCounterIndex=0
#  dGlassDimX=2.            # [cm]
#  dGlassDimY=2.            # [cm]
#  iModuleType=9
#  iModuleIndex=0
#  iCounterIndex=0
#  dGlassDimX=25.           # [cm]
#  dGlassDimY=26.7          # [cm]
  iModuleType=7
  iModuleIndex=0
  iCounterIndex=0
  dGlassDimX=32.            # [cm]
  dGlassDimY=27.            # [cm]
  dHitDistance=0.01        # [cm] avalanche impact parameter
  dRelaxationTime=5.0      # [s] MRPC relaxation time constant
  dWorkingCoefficient=1.5  # charge distribution scale factor
  bPlaneVertex=0           # NOTE: takes no effect in this run mode
  bEquidistantEvents=0     # 1: event time intervals are fixed
                           # 0: event time intervals are exponentially distributed

  bRunReco=1
  bRunQA=1
  bAlternativeBranchNames=0
  bBuildCbmEvents=0
  bBuildHits=1
  bHitCompatibility=0
  bSpotResponse=1
  bNoTransportPhysics=1
  bTwoSidedReadout=1
  bCalibrateDigis=1
  bStorePointsInTS=1

fi
################################################################################
# HAF-specific settings                                                        #
################################################################################
if [ "${cRunMode}" == "haf" ]; then
  bSingleCounter=0
  dParticleFlux=1000.      # [Hz/cm^2]
  iModuleType=7
  iModuleIndex=0
  iCounterIndex=0
  dGlassDimX=32.           # [cm]
  dGlassDimY=27.           # [cm]
  dHitDistance=0.          # NOTE: takes no effect in this run mode
  dRelaxationTime=0.       # NOTE: takes no effect in this run mode
  dWorkingCoefficient=0.   # NOTE: takes no effect in this run mode
  bPlaneVertex=0
  bEquidistantEvents=0     # 1: event time intervals are fixed
                           # 0: event time intervals are exponentially distributed

  iProductionMode=0        # 0: calibrated digi and hit output (with reco, with QA)
                           # 1: uncalibrated digi output (no reco, no QA)
                           # 2: calibrated digi output assigned to 'CbmEvent' objects (no hits, no QA)
                           # 3: uncalibrated digi output assigned to 'CbmEvent' objects (no hits, no QA)

  bHitCompatibility=1
  bSpotResponse=0
  bNoTransportPhysics=0

################################################################################
# HAF-specific production settings                                             #
################################################################################
  if [ "${iProductionMode}" == "0" ]; then
    bRunReco=1
    bRunQA=1
    bAlternativeBranchNames=1
    bBuildCbmEvents=0
    bBuildHits=1
    bTwoSidedReadout=1
    bCalibrateDigis=1
    bStorePointsInTS=1
  elif [ "${iProductionMode}" == "1" ]; then
    bRunReco=0
    bRunQA=0
    bAlternativeBranchNames=1
    bBuildCbmEvents=0
    bBuildHits=0
    bTwoSidedReadout=0
    bCalibrateDigis=0
    bStorePointsInTS=1
  elif [ "${iProductionMode}" == "2" ]; then
    bRunReco=1
    bRunQA=0
    bAlternativeBranchNames=0
    bBuildCbmEvents=1
    bBuildHits=0
    bTwoSidedReadout=0
    bCalibrateDigis=1
    bStorePointsInTS=0
  elif [ "${iProductionMode}" == "3" ]; then
    bRunReco=1
    bRunQA=0
    bAlternativeBranchNames=0
    bBuildCbmEvents=1
    bBuildHits=0
    bTwoSidedReadout=0
    bCalibrateDigis=0
    bStorePointsInTS=0
  else
    echo "HAF production mode ${iProductionMode} not implemented"
    return 1
  fi
fi

################################################################################
# HIC-specific general settings                                                #
################################################################################
if [ "${cRunMode}" == "hic" ]; then
  cSystem=arpb
  dBeamMomentum=19
  cCentrality=mbias
  dBeamIntensity=17.6e6    # [Hz]
  dTargetInteraction=0.04  # range [0.,1.) here: 4% (equivalent to 2 mm of Pb)
                           #                     8% (equivalent to 4 mm of Pb)
                           # IMPORTANT NOTE: The product of the preceding two
                           #                 parameters corresponds to the
                           #                 event rate in the setup.
  bGenerateDarkRate=1      # 0: do not generate any counter dark rate (ignore all 'CounterDarkRate'
                           #    fields in parameters/tof/tof_${cSetup}.digitb.par)
  bIgnoreInterference=0    # 1: treat overlapping digital signals independently
                           # 0: merge overlapping digital signals
  bIgnoreDeadTime=0        # 1: assume dead-time free TDCs (ignore all 'ChannelDeadTime' fields
                           #    in parameters/tof/tof_${cSetup}.digitb.par)
  iMCTimeLimit=0-20:00:00  # expected MC job completion time (TIMEOUT if exceeded)
  iNEventsPerFile=10000
  iNParallelMCRuns=10      # WARNING: Make sure that this number divides without rest
                           #          the required numbers of MC events/files per
                           #          digitization run!
  bReuseMCFiles=1

  iProductionMode=0        # 0: calibrated digi and hit output (with reco, with QA)
                           # 1: uncalibrated digi output (no reco, no QA)
                           # 2: calibrated digi output assigned to 'CbmEvent' objects (no hits, no QA)
                           # 3: uncalibrated digi output assigned to 'CbmEvent' objects (no hits, no QA)

  bHitCompatibility=1
  bSpotResponse=0
  bNoTransportPhysics=0

################################################################################
# HIC-specific production settings                                             #
################################################################################
  if [ "${iProductionMode}" == "0" ]; then
    bRunReco=1
    bRunQA=1
    bAlternativeBranchNames=1
    bBuildCbmEvents=0
    bBuildHits=1
    bTwoSidedReadout=1
    bCalibrateDigis=1
    bStorePointsInTS=1
  elif [ "${iProductionMode}" == "1" ]; then
    bRunReco=0
    bRunQA=0
    bAlternativeBranchNames=1
    bBuildCbmEvents=0
    bBuildHits=0
    bTwoSidedReadout=0
    bCalibrateDigis=0
    bStorePointsInTS=1
  elif [ "${iProductionMode}" == "2" ]; then
    bRunReco=1
    bRunQA=0
    bAlternativeBranchNames=0
    bBuildCbmEvents=1
    bBuildHits=0
    bTwoSidedReadout=0
    bCalibrateDigis=1
    bStorePointsInTS=0
  elif [ "${iProductionMode}" == "3" ]; then
    bRunReco=1
    bRunQA=0
    bAlternativeBranchNames=0
    bBuildCbmEvents=1
    bBuildHits=0
    bTwoSidedReadout=0
    bCalibrateDigis=0
    bStorePointsInTS=0
  else
    echo "HIC production mode ${iProductionMode} not implemented"
    return 1
  fi
fi

################################################################################
# Debug settings (overwrite parts of the above)                                #
################################################################################
#iNLogicalCores=7
#iTimeLimit=0-01:00:00
#iNParallelSpills=5
#dBeamMomentum=30
#dBeamIntensity=1.1e4
#iNParallelMCRuns=3
#iMCTimeLimit=0-00:30:00
#bTwoSidedReadout=1
#bRunWithMemory=0
#bGenerateDarkRate=0
#bIgnoreInterference=1
#bIgnoreDeadTime=1
#bReuseMCFiles=0
#bStorePointsInTS=0
#bRunReco=1
#bRunQA=0
#bHitCompatibility=1
#bNoTransportPhysics=1
################################################################################
################################################################################

macroDir=/lustre/nyx/cbm/users/csimon/cbmroot/development/csimon/src/macro/tof/digitizer
#mcDir=/lustre/nyx/cbm/prod/gen/urqmd/${cSystem}/${dBeamMomentum}gev/${cCentrality}
mcDir=/lustre/nyx/cbm/users/csimon/prod/gen/urqmd/${cSystem}/${dBeamMomentum}gev/${cCentrality}
mcFileName=${mcDir}/urqmd.${cSystem}.${dBeamMomentum}gev.${cCentrality}

wdir=/lustre/nyx/cbm/users/csimon/simulation/production/${cSetup}


if [ "${cRunMode}" == "scf" ]; then

  if [ ! -z "$1" ]; then
    dParticleFlux=$1
  fi

  if [ ! -z "$2" ]; then
    dRelaxationTime=$2
  fi

  if [ ! -z "$3" ]; then
    dWorkingCoefficient=$3
  fi

  if [ ! -z "$4" ]; then
    dHitDistance=$4
  fi

  if [ ! -z "$5" ]; then
    iNLogicalCores=$5
  fi

  if [ ! -z "$6" ]; then
    iTimeLimit=$6
  fi


  if [ "${iNSerialSpills}" == "1" ]; then
    dSpillBreakLength=0.
    bAverageSerialSpills=0
  fi

  if [ "${bAverageSerialSpills}" == "1" ]; then
  # if serials spills are supposed to be averaged over, distinguishing between different
  # spill break lengths is redundant
    dSpillBreakLength=5.
  fi

  wdir=${wdir}/rm_scf
  wdir=${wdir}/${iModuleType}_${iModuleIndex}_${iCounterIndex}

  wdir=${wdir}/`echo "(${dRelaxationTime}*1000)/1" | bc`ms_relaxation
  wdir=${wdir}/`echo "(${dWorkingCoefficient}*100)/1" | bc`pct_charge
  wdir=${wdir}/`echo "(${dHitDistance}*10000)/1" | bc`um_avalanche # FIXME: us_avalanche

  if [ "${bSpotResponse}" == "0" ]; then
    wdir=${wdir}/full
  else
    wdir=${wdir}/spot
  fi

  wdir=${wdir}/${dParticleFlux%?}Hz_per_sq_cm
  wdir=${wdir}/${dSpillLength%?}s_spill
  wdir=${wdir}/${iNSerialSpills}spills_per_file

  if [ "${bAverageSerialSpills}" == "1" ]; then
    wdir=${wdir}/average
  else
    wdir=${wdir}/ongoing
  fi

  wdir=${wdir}/${dSpillBreakLength%?}s_break

elif [ "${cRunMode}" == "haf" ]; then

  wdir=${wdir}/rm_haf
  wdir=${wdir}/${iModuleType}_${iModuleIndex}_${iCounterIndex}

  if [ "${bSpotResponse}" == "0" ]; then
    wdir=${wdir}/full
  else
    wdir=${wdir}/spot
  fi

  if [ "${bPlaneVertex}" == "0" ]; then
    wdir=${wdir}/origin
  else
    wdir=${wdir}/plane
  fi

  wdir=${wdir}/${dParticleFlux%?}Hz_per_sq_cm

  if [ "${bRunWithMemory}" == "1" ]; then
    wdir=${wdir}/memory
  else
    wdir=${wdir}/normal
  fi

  wdir=${wdir}/${dSpillLength%?}s_spill
  wdir=${wdir}/${dSpillBreakLength%?}s_break
  wdir=${wdir}/${iNSerialSpills}spills_per_file

  if [ "${bAverageSerialSpills}" == "1" ]; then
    wdir=${wdir}/average
  else
    wdir=${wdir}/ongoing
  fi

  wdir=${wdir}/pm_${iProductionMode}

  wdir=${wdir}/`date -u +"%Y-%m-%dT%H-%M-%SZ"`

elif [ "${cRunMode}" == "hic" ]; then

  if [ ! -d "${mcDir}" ]; then
    echo 'MC input directory does not exist'
    return 1
  fi

  # bc cannot handle scientific e-notation
  dBeamIntensityBC=`echo ${dBeamIntensity} | sed -e 's/[eE]+*/\\*10\\^/'`

  dEventRate=`echo "${dBeamIntensityBC}*${dTargetInteraction}" | bc -l`
  dSpillTime=`echo "${dSpillLength}*${iNSerialSpills}" | bc -l`

  # number of MC events needed for a single digitization run
  iNEvents=`echo "1/2*(2*${dEventRate}*${dSpillTime} +5*sqrt(4*${dEventRate}*${dSpillTime}+25) +25) +0.5" | bc -l`
  iNEvents=`echo "${iNEvents}/1" | bc`
  # Time intervals between events are exponentially distributed. The period of time
  # covered by n subsequent events is described by an Erlang distribution of shape n.
  # The above calculation results in the required minimum number of events to be simulated for
  # a given number of serial spills. With a safety margin of 5 Erlang standard deviations
  # the hypothetical next spill in line which is not meant to be processed is reached
  # by 'CbmFileSource'. If the amount of events were insufficient to cross the final
  # spill break, no dark-rate digis would be produced therein.

  iNMCArrayJobs=$((iNParallelSpills*iNParallelMCRuns))
  iNMCFiles=$(((iNEvents+iNEventsPerFile-1)/iNEventsPerFile))
  iNFilesPerMCRun=$(((iNMCFiles+iNParallelMCRuns-1)/iNParallelMCRuns))

  iNAvailMCFiles=`find ${mcDir} -maxdepth 1 -name '*[0-9][0-9][0-9][0-9][0-9].root' | wc -l`

  if [ "${bReuseMCFiles}" == "0" ]; then
    iNMCFiles=$((iNMCFiles*iNParallelSpills))
  fi

  if [ "${iNAvailMCFiles}" -lt "${iNMCFiles}" ]; then
    echo 'not enough MC input files available'
    return 1
  fi

  # check if enough MC files are contiguously indexed
  MCFileList=( $(find ${mcDir} -maxdepth 1 -name '*[0-9][0-9][0-9][0-9][0-9].root' -printf '%f\n' | rev | cut -c6-10 | rev | sort -n) )
  MCFileList=("${MCFileList[@]:0:iNMCFiles}")
  iNMissingMCFiles=`join -v1 <(seq -w ${MCFileList[-1]}) <(printf '%s\n' ${MCFileList[@]}) | wc -l`

  if [ "${iNMissingMCFiles}" != "0" ]; then
    echo 'MC input files are not contiguously indexed'
    echo "MC file indices missing: ${iNMissingMCFiles}"
    join -v1 <(seq -w ${MCFileList[-1]}) <(printf '%s\n' ${MCFileList[@]})
    return 1
  fi

  iNEventsPerMCRun=$(((iNEvents+iNParallelMCRuns-1)/iNParallelMCRuns))

  wdir=${wdir}/rm_hic
  wdir=${wdir}/${cSystem}
  wdir=${wdir}/${dBeamMomentum}gev
  wdir=${wdir}/${cCentrality}

  if [ "${bRunWithMemory}" == "1" ]; then
    wdir=${wdir}/memory
  else
    wdir=${wdir}/normal
  fi

  wdir=${wdir}/${dSpillLength%?}s_spill
  wdir=${wdir}/${dSpillBreakLength%?}s_break
  wdir=${wdir}/${iNSerialSpills}spills_per_file

  if [ "${bAverageSerialSpills}" == "1" ]; then
    wdir=${wdir}/average
  else
    wdir=${wdir}/ongoing
  fi

  wdir=${wdir}/pm_${iProductionMode}

  wdir=${wdir}/`date -u +"%Y-%m-%dT%H-%M-%SZ"`

fi


mkdir -p ${wdir}

if [ ! -d "${wdir}" ]; then
  echo 'working directory does not exist and cannot be created'
  return 1
fi


iNPhysicalCores=`echo "${iNLogicalCores}/2 + 1" | bc`

if [ "${bRunWithMemory}" == "0" ]; then
  iNPhysicalCores=1
fi


if [ "${cRunMode}" == "scf" -o "${cRunMode}" == "haf" ]; then

  cDep=$(sbatch --time=${iTimeLimit} --cpus-per-task=${iNPhysicalCores} --partition=${cPartition} --workdir=${wdir} --array=1-${iNParallelSpills} ${macroDir}/run_flood.sh \
         ${cSetup} ${iNSerialSpills} ${dSpillLength} ${dSpillBreakLength} ${dParticleFlux} ${iModuleType} ${iModuleIndex} ${iCounterIndex} \
         ${dGlassDimX} ${dGlassDimY} ${bSpotResponse} ${bAverageSerialSpills} ${dHitDistance} ${dRelaxationTime} ${bSingleTimeSlice} \
         ${bSingleCounter} ${bPlaneVertex} ${bHitCompatibility} ${bEquidistantEvents} ${bRunWithMemory} ${bRunQA} ${cGeometry} ${bAlternativeBranchNames} \
         ${bBuildCbmEvents} ${bBuildHits} ${bRunReco} ${bNoTransportPhysics} ${bTwoSidedReadout} ${bCalibrateDigis} ${bStorePointsInTS} ${dWorkingCoefficient})
  cDep=${cDep##* }

elif [ "${cRunMode}" == "hic" ]; then

  cDepMC=$(sbatch --time=${iMCTimeLimit} --partition=${cPartition} --workdir=${wdir} --array=1-${iNMCArrayJobs} ${macroDir}/run_hic_mc.sh \
           ${cSetup} ${bNoTransportPhysics} ${iNEventsPerMCRun} ${iNFilesPerMCRun} ${iNParallelMCRuns} ${bReuseMCFiles} \
           ${mcFileName} ${cGeometry} ${dTargetInteraction})
  cDepMC=${cDepMC##* }

  cDep=$(sbatch --dependency=afterany:${cDepMC} --time=${iTimeLimit} --cpus-per-task=${iNPhysicalCores} --partition=${cPartition} --workdir=${wdir} --array=1-${iNParallelSpills} \
         ${macroDir}/run_hic.sh \
         ${cSetup} ${iNSerialSpills} ${dSpillLength} ${dSpillBreakLength} ${dBeamIntensity} \
         ${bAverageSerialSpills} ${bSingleTimeSlice} ${bCalibrateDigis} ${iNEvents} ${iNParallelMCRuns} ${cGeometry} ${dTargetInteraction} ${bTwoSidedReadout} \
         ${bRunWithMemory} ${bGenerateDarkRate} ${bIgnoreInterference} ${bIgnoreDeadTime} ${bStorePointsInTS} ${bRunQA} ${bRunReco} ${bHitCompatibility} \
         ${bAlternativeBranchNames} ${bBuildCbmEvents} ${bBuildHits})
  cDep=${cDep##* }

else

    echo "run mode ${cRunMode} not implemented"
    return 1

fi


if [ "${bRunQA}" == "1" ]; then

  sbatch --dependency=afterany:${cDep} --partition=${cPartition} --workdir=${wdir} ${macroDir}/combine_histograms.sh ${iNParallelSpills} ${cGeometry} ${bSpotResponse}

fi