diff --git a/MQ/hitbuilder/startMQ_Mcbm2018.sh b/MQ/hitbuilder/startMQ_Mcbm2018.sh
new file mode 100644
index 0000000000000000000000000000000000000000..9232587e7c3109ac0571dc95e105a1bec5608031
--- /dev/null
+++ b/MQ/hitbuilder/startMQ_Mcbm2018.sh
@@ -0,0 +1,286 @@
+#!/bin/bash
+$FAIRROOTPATH/bin/shmmonitor --cleanup
+
+if [ -z "$1" ]; then 
+    _runname="359"
+else 
+    _runname=$1
+fi
+ 
+if [ -z "$2" ]; then 
+_reqmod=-193
+else 
+_reqmod=$2
+fi
+
+if [ -z "$3" ]; then 
+_pulmode=1
+else 
+_pulmode=$3
+fi
+
+if [ -z "$4" ]; then 
+    _reqtint=100
+else
+    _reqtint=$4 
+fi
+
+_Opt=$5
+
+_ntimeslices=-1
+#_ntimeslices=10000
+_iUnp=4
+_batch=1
+_pulmulmin=5
+_pultotmin=50
+_pultotmax=500
+_puldetref=12 # TSR=022
+#_puldetref=16 # TSR=032
+#_puldetref=17 # TSR=032
+
+#_tofftof=0.  
+_tofftof=30. 
+
+if [[ $_reqmod -eq 1 ]]; then
+  _iUnp=1
+fi
+
+#_dirname=/lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/input/$_runname/
+_outdir=/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/data/
+#_datapath=~/KRONOS/CBM/cbmroot/trunk
+#_dirname=$_datapath/macro/beamtime/mcbm2018/input/$_runname/
+#_outdir=$_datapath/macro/beamtime/mcbm2018/data/
+ 
+#_filename="./tsaData/${_runname}_pn05_*.tsa;./tsaData/${_runname}_pn06_*.tsa;./tsaData/${_runname}_pn07_*.tsa"
+_filename="./tsaData/${_runname}_pn02_*.tsa;./tsaData/${_runname}_pn04_*.tsa;./tsaData/${_runname}_pn05_*.tsa;./tsaData/${_runname}_pn06_*.tsa;./tsaData/${_runname}_pn07_*.tsa;./tsaData/${_runname}_pn08_*.tsa;./tsaData/${_runname}_pn10_*.tsa"
+#_filename="./tsaData2020/${_runname}_pn02_*.tsa;./tsaData2020/${_runname}_pn04_*.tsa;./tsaData2020/${_runname}_pn05_*.tsa;./tsaData2020/${_runname}_pn06_*.tsa;./tsaData2020/${_runname}_pn08_*.tsa;./tsaData2020/${_runname}_pn10_*.tsa;./tsaData2020/${_runname}_pn11_*.tsa;./tsaData2020/${_runname}_pn12_*.tsa;./tsaData2020/${_runname}_pn13_*.tsa;./tsaData2020/${_runname}_pn15_*.tsa"
+_digifile=$_runname.$_reqtint.$_reqmod.${_pulmode}${_Opt}.root
+
+# ASCII files  
+#_mapfile=/lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/MQ/hitbuilder/MapTofGbtx.par
+_mapfile=/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/mTofPar.par
+if [ $_runname -ge 707 ] && [ $_runname -lt 754 ]; then 
+  _mapfile=/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2020/mTofPar_2Stack.par
+else
+  if [ $_runname -ge 754 ]; then 
+    _mapfile=/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2020/mTofPar_3Stack.par
+    _puldetref=12 # TSR=022
+  fi
+fi 
+#_mapfile=/lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/etofPar.par
+#_digibdffile=/lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/parameters/tof/v18j_cosmicHD.digibdf.par
+#_digiparfile=/lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/parameters/tof/tof_v18j_cosmicHD.digi.par
+_digibdffile=/lustre/cbm/users/nh/CBM/cbmroot/trunk/parameters/tof/v19b_mcbm.digibdf.par
+_digiparfile=/lustre/cbm/users/nh/CBM/cbmroot/trunk/parameters/tof/tof_v19b_mcbm.digi.par
+
+# ROOT files 
+#_geofile=/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/tof_v18l_mCbm.par.root
+_geofile=/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/tof_mcbm_beam_2019_03.par.root
+
+# MQ ports
+_pPar=5603
+_pSam=5613
+_pCmd=5623
+_pDig=5633
+
+rm -v nohup.out 
+rm -v *log
+rm all_*
+
+PARAMETERSERVER="parmq-server"
+echo pkill $PARAMETERSERVER
+pkill -9 $PARAMETERSERVER
+pkill -9 $PARAMETERSERVER
+sleep 1
+PARAMETERSERVER+=" --id parmq-server"
+PARAMETERSERVER+=" --channel-name parameters"
+PARAMETERSERVER+=" --channel-config name=parameters,type=rep,method=bind,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:$_pPar"
+#PARAMETERSERVER+=" --libs-to-load libCbmTof;libCbmFlibMcbm2018"
+PARAMETERSERVER+=" --first-input-name $_mapfile;$_digiparfile;$_digibdffile"
+PARAMETERSERVER+=" --first-input-type ASCII"
+PARAMETERSERVER+=" --second-input-name $_geofile"
+PARAMETERSERVER+=" --second-input-type ROOT"
+PARAMETERSERVER+=" --severity INFO"
+if  [[ $_batch = 1 ]]; then 
+PARAMETERSERVER+=" --control static"
+PARAMETERSERVER+=" --log-to-file ParServ.out"
+nohup /lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/build/bin/MQ/parmq/$PARAMETERSERVER &
+else
+xterm -geometry 80x23+0+340 -hold -e /lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/build/bin/MQ/parmq/$PARAMETERSERVER &
+#xterm -geometry 80x23+500+0 -hold -e /cvmfs/fairroot.gsi.de/fairroot/v18.0.6_fairsoft-may18/bin/$PARAMETERSERVER &
+fi
+
+SAMPLER="TsaMultiSamplerTof"
+SAMPLER+=" --id sampler1"
+SAMPLER+=" --max-timeslices $_ntimeslices"
+#SAMPLER+=" --max-timeslices 1000000"
+#SAMPLER+=" --flib-port 10"
+#SAMPLER+=" --dirname $_dirname"
+SAMPLER+=" --filename $_filename"
+#SAMPLER+=" --flib-host myHost"
+SAMPLER+=" --channel-config name=tofcomponent,type=push,method=bind,rateLogging=0,transport=zeromq,address=tcp://*:$_pSam"
+SAMPLER+=" --channel-config name=syscmd,type=pub,method=bind,rateLogging=0,transport=zeromq,address=tcp://*:$_pCmd"
+#SAMPLER+=" --transport shmem"
+#SAMPLER+=" --transport zeromq"
+#SAMPLER+=" --transport nanomsg"
+#SAMPLER+=" --severity WARN"
+SAMPLER+=" --severity INFO"
+#SAMPLER+=" --severity DEBUG"
+SAMPLER+=" --SelectComponents 1"
+if  [[ $_batch = 1 ]]; then 
+SAMPLER+=" --log-to-file Sampl.out"
+SAMPLER+=" --control static"
+nohup  /lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/build/bin/MQ/source/$SAMPLER & 
+else 
+xterm -geometry 80x23+0+0 -hold -e /lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/build/bin/MQ/source/$SAMPLER &
+fi
+ 
+while (( _iUnp > 0 )); do
+
+UNPACKER="UnpackTofMcbm2018"
+UNPACKER+=" --id unpack$_iUnp"
+UNPACKER+=" --channel-config name=tofcomponent,type=pull,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:$_pSam"
+UNPACKER+=" --channel-config name=parameters,type=req,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:$_pPar"
+UNPACKER+=" --channel-config name=tofdigis,type=push,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:$_pDig"
+UNPACKER+=" --channel-config name=syscmd,type=sub,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:$_pCmd"
+#UNPACKER+=" --transport shmem"
+#UNPACKER+=" --severity DEBUG"
+UNPACKER+=" --severity  INFO"
+UNPACKER+=" --SelectComponents 1"
+#UNPACKER+=" --ReqBeam      20486" # diamond -> 0x00005006 v14a
+UNPACKER+=" --ReqBeam      10246" # diamond -> 0x00002806 v21a
+if  [[ $_reqmod -lt 1 ]]; then
+    UNPACKER+=" --ReqMode 0"
+    case $_reqmod in
+	0)
+	    ;;
+       -1)
+	    UNPACKER+=" --ReqDet0       20486" # diamond -> 0x00005006
+	    UNPACKER+=" --ReqDet1       65590" # RPC 031 -> 0x00010036
+	    UNPACKER+=" --ReqDet2       65606" # RPC 041
+	    ;;
+       -2) 
+	    UNPACKER+=" --ReqDet0       20486" # diamond
+	    UNPACKER+=" --ReqDet1      131126" # RPC 032
+	    UNPACKER+=" --ReqDet2      131142" # RPC 042
+	    ;; # for ceramics
+
+       -3)
+	    UNPACKER+=" --ReqDet0       20486" # diamond
+	    UNPACKER+=" --ReqDet1      196662" # RPC 033
+	    UNPACKER+=" --ReqDet2      196678" # RPC 043
+	    ;; # for BUC
+	    
+       -4) # v21a address mode
+	    UNPACKER+=" --ReqDet0       10246" # diamond
+	    UNPACKER+=" --ReqDet1       65542" # RPC 022       
+	    UNPACKER+=" --ReqDet2       65574" # RPC 022       
+	    ;;
+       
+       -190) 
+	    UNPACKER+=" --ReqDet0       20486"  # diamond
+	    UNPACKER+=" --ReqDet1       65606"  # RPC 041
+	    UNPACKER+=" --ReqDet2       36870"  # RPC 900
+	    UNPACKER+=" --ReqDet3      102406"  # RPC 901
+	    ;; # for  double stack calibration
+       
+       -191) 
+	    UNPACKER+=" --ReqDet0       20486" # diamond
+	    UNPACKER+=" --ReqDet1       65606"  # RPC 041
+	    UNPACKER+=" --ReqDet2       24582"  # RPC 600
+	    UNPACKER+=" --ReqDet3       90118"   # RPC 601
+	    ;; # for USTC counter analysis
+
+       -192) 
+	    UNPACKER+=" --ReqDet0       65606"  # RPC 041
+	    UNPACKER+=" --ReqDet1       36870"  # RPC 900
+	    UNPACKER+=" --ReqDet2     102406"   # RPC 901
+	    ;; # for  double stack calibration
+
+       -193) BeamBeam
+	    UNPACKER+=" --ReqDet0       65606"  # RPC 041
+	    UNPACKER+=" --ReqDet1       24582"  # RPC 600
+	    UNPACKER+=" --ReqDet2       90118"  # RPC 601
+	    ;; # for USTC counter analysis
+
+       -194) 
+	    UNPACKER+=" --ReqDet0       20486"  # diamond
+	    UNPACKER+=" --ReqDet1       24582"  # RPC 600
+	    UNPACKER+=" --ReqDet2       90118"  # RPC 601
+	    ;; # for USTC counter analysis
+	    
+       *)
+	    echo ReqMode $_reqmod not yet defined, exiting
+	    exit 1
+	    ;;
+    esac;
+else
+    UNPACKER+=" --ReqMode $_reqmod"
+fi
+UNPACKER+=" --ReqTint $_reqtint"
+UNPACKER+=" --PulserMode $_pulmode"
+UNPACKER+=" --PulMulMin $_pulmulmin"
+UNPACKER+=" --PulTotMin $_pultotmin"
+UNPACKER+=" --PulTotMax $_pultotmax"
+UNPACKER+=" --ToffTof $_tofftof"
+UNPACKER+=" --RefModType    5"
+UNPACKER+=" --RefModId      0"
+UNPACKER+=" --RefCtrType    0"
+UNPACKER+=" --RefCtrId      0"
+if  [[ $_batch = 1 ]]; then 
+UNPACKER+=" --control static"
+UNPACKER+=" --log-to-file Unp$_iUnp.out"
+#echo nohup $UNPACKER 
+nohup /lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/build/bin/MQ/unpacker/$UNPACKER &
+else 
+xterm -geometry 110x23+520+0 -hold -e /lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/build/bin/MQ/unpacker/$UNPACKER & 
+fi
+
+(( _iUnp -= 1 ))
+done 
+
+HITBUILDER="HitBuilderTof"
+HITBUILDER+=" --id hitbuilder1"
+HITBUILDER+=" --channel-config name=tofdigis,type=pull,method=bind,rateLogging=0,transport=zeromq,address=tcp://*:$_pDig"
+HITBUILDER+=" --channel-config name=parameters,type=req,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:$_pPar"
+HITBUILDER+=" --channel-config name=syscmd,type=sub,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:$_pCmd"
+#HITBUILDER+=" --channel-config name=tofhits,type=push,method=bind,transport=shmem,address=tcp://127.0.0.1:5557"
+#HITBUILDER+=" --channel-config name=tofcalib,type=push,method=bind,transport=shmem,address=tcp://127.0.0.1:5558"
+#HITBUILDER+=" --transport shmem"
+#HITBUILDER+=" --severity DEBUG"
+HITBUILDER+=" --severity INFO"
+#HITBUILDER+=" --severity WARN"
+HITBUILDER+=" --OutRootFile $_outdir$_digifile"
+#HITBUILDER+=" --MaxEvent 10000000"
+#HITBUILDER+=" --RunId 1552883952"
+HITBUILDER+=" --RunId 1601311083"
+HITBUILDER+=" --PulserMode $_pulmode"
+HITBUILDER+=" --PulMulMin $_pulmulmin"
+HITBUILDER+=" --PulTotMin $_pultotmin"
+HITBUILDER+=" --PulTotMax $_pultotmax"
+HITBUILDER+=" --PulDetRef $_puldetref"
+HITBUILDER+=" --ReqTint $_reqtint"
+#HITBUILDER+=" --ReqBeam      20486" # diamond -> 0x00005006
+HITBUILDER+=" --BRefType  5"
+HITBUILDER+=" --BRefSm    0"
+HITBUILDER+=" --BRefDet   0"
+HITBUILDER+=" --DutType   0"
+HITBUILDER+=" --DutSm     3"
+HITBUILDER+=" --DutRpc    1"
+HITBUILDER+=" --SelType   0"
+HITBUILDER+=" --SelSm     4"
+HITBUILDER+=" --SelRpc    1"
+HITBUILDER+=" --Sel2Type  5"
+HITBUILDER+=" --Sel2Sm    0"
+HITBUILDER+=" --Sel2Rpc   0"
+if [[ $_reqmod -eq 1 ]]; then
+HITBUILDER+=" --Mode   1"
+fi
+
+if  [[ $_batch = 1 ]]; then 
+HITBUILDER+=" --control static"
+HITBUILDER+=" --log-to-file HitBuild.out"
+nohup  /lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/build/bin/MQ/hitbuilder/$HITBUILDER &
+else
+xterm -geometry 120x23+1400+0 -hold -e /lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/build/bin/MQ/hitbuilder/$HITBUILDER &
+fi
diff --git a/MQ/hitbuilder/startMQ_cosmicHD.sh b/MQ/hitbuilder/startMQ_cosmicHD.sh
new file mode 100755
index 0000000000000000000000000000000000000000..1876f84799cf0efc0b16c235dfc4182ab1acca7a
--- /dev/null
+++ b/MQ/hitbuilder/startMQ_cosmicHD.sh
@@ -0,0 +1,176 @@
+#!/bin/bash
+$FAIRROOTPATH/bin/shmmonitor --cleanup
+
+if [ -z "$1" ]; then 
+    _runname="HD_cosmic_`date +%F`_`date +%T`"
+else 
+    _runname=$1
+fi
+
+_reqmod=3
+_reqtint=50
+_ntimeslices=-1
+_iUnp=4
+_batch=0
+_pulmode=0
+_pulmulmin=3
+_pultotmin=50
+_pultotmax=500
+_puldetref=16 # TSR=032
+#_puldetref=17 # TSR=032
+
+#_tofftof=0.  
+_tofftof=-30. 
+_version=0
+
+if [[ $_reqmod -eq 1 ]]; then
+  _iUnp=1
+fi
+
+_dirname=/home/cbm/software/cbmroot/git/cbmroot/macro/beamtime/hd2020/input/$_runname/
+#_outdir=/home/cbm/software/cbmroot/git/cbmroot/macro/beamtime/hd2020/data/
+_outdir=/d/cbm/HD2020/
+#_datapath=~/KRONOS/CBM/cbmroot/git/cbmroot
+_filename=$_runname*.tsa
+#_filename=/home/cbm/software/cbmroot/git/cbmroot/macro/beamtime/hd2020/input/$_runname/*.tsa
+#_outdir=$_datapath/macro/beamtime/mcbm2018/data/
+ 
+#_filename="./tsaData/_pn05_*.tsa;./tsaData/_pn06_*.tsa;./tsaData/_pn07_*.tsa"
+_digifile=$_runname.$_reqtint.$_reqmod.$_pulmode.$_version.root
+
+# ASCII files  
+_mapfile=/home/cbm/software/cbmroot/git/cbmroot/macro/beamtime/hd2020/mTofParBuc.par  
+_digibdffile=/home/cbm/software/cbmroot/git/cbmroot/parameters/tof/tof_v20b_cosmicHD.digibdf.par
+_digiparfile=/home/cbm/software/cbmroot/git/cbmroot/parameters/tof/tof_v20b_cosmicHD.digi.par
+
+# ROOT files 
+_geofile=/home/cbm/software/cbmroot/git/cbmroot/macro/beamtime/hd2020/tof_cosmicHD_2020_01.par.root
+
+rm -v nohup.out 
+rm -v *log
+rm all_*
+
+PARAMETERSERVER="parmq-server"
+echo pkill $PARAMETERSERVER
+pkill -9 $PARAMETERSERVER
+sleep 1
+PARAMETERSERVER+=" --id parmq-server"
+PARAMETERSERVER+=" --channel-name parameters"
+PARAMETERSERVER+=" --channel-config name=parameters,type=rep,method=bind,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:5005"
+#PARAMETERSERVER+=" --libs-to-load libCbmTof;libCbmFlibMcbm2018"
+PARAMETERSERVER+=" --first-input-name $_mapfile;$_digiparfile;$_digibdffile"
+#PARAMETERSERVER+=" --first-input-name $_mapfile;$_digibdffile"
+PARAMETERSERVER+=" --first-input-type ASCII"
+PARAMETERSERVER+=" --second-input-name $_geofile"
+PARAMETERSERVER+=" --second-input-type ROOT"
+PARAMETERSERVER+=" --severity INFO"
+if  [[ $_batch = 1 ]]; then 
+PARAMETERSERVER+=" --control static"
+PARAMETERSERVER+=" --log-to-file ParServ.out"
+nohup /home/cbm/software/cbmroot/git/cbmroot/build/bin/MQ/parmq/$PARAMETERSERVER &
+else
+xterm -geometry 80x23+0+340 -hold -e /home/cbm/software/cbmroot/git/cbmroot/build/bin/MQ/parmq/$PARAMETERSERVER &
+#xterm -geometry 80x23+500+0 -hold -e /home/cbm/starsoft/fairroot_v18.0.6-fairsoft_may18p1_root6/bin/$PARAMETERSERVER &
+fi
+
+SAMPLER="TsaMultiSamplerTof"
+SAMPLER+=" --id sampler1"
+SAMPLER+=" --max-timeslices $_ntimeslices"
+#SAMPLER+=" --max-timeslices 1000000"
+SAMPLER+=" --flib-port 5556"
+SAMPLER+=" --flib-host cbmin008"
+#SAMPLER+=" --dirname $_dirname"
+#SAMPLER+=" --filename $_filename"
+SAMPLER+=" --channel-config name=tofcomponent,type=push,method=bind,rateLogging=0,transport=zeromq,address=tcp://*:5655"
+SAMPLER+=" --channel-config name=syscmd,type=pub,method=bind,rateLogging=0,transport=zeromq,address=tcp://*:5666"
+#SAMPLER+=" --transport shmem"
+#SAMPLER+=" --transport zeromq"
+#SAMPLER+=" --transport nanomsg"
+#SAMPLER+=" --severity WARN"
+SAMPLER+=" --severity INFO"
+#SAMPLER+=" --severity DEBUG"
+SAMPLER+=" --SelectComponents 0"
+if  [[ $_batch = 1 ]]; then 
+SAMPLER+=" --log-to-file Sampl.out"
+SAMPLER+=" --control static"
+nohup  /home/cbm/software/cbmroot/git/cbmroot/build/bin/MQ/source/$SAMPLER & 
+else 
+xterm -geometry 80x23+0+0 -hold -e /home/cbm/software/cbmroot/git/cbmroot/build/bin/MQ/source/$SAMPLER &
+fi
+ 
+while (( _iUnp > 0 )); do
+
+UNPACKER="UnpackTofMcbm2018"
+UNPACKER+=" --id unpack$_iUnp"
+UNPACKER+=" --channel-config name=tofcomponent,type=pull,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:5655"
+UNPACKER+=" --channel-config name=parameters,type=req,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:5005"
+UNPACKER+=" --channel-config name=tofdigis,type=push,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:5656"
+UNPACKER+=" --channel-config name=syscmd,type=sub,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:5666"
+#UNPACKER+=" --transport shmem"
+#UNPACKER+=" --severity DEBUG"
+UNPACKER+=" --severity  INFO"
+UNPACKER+=" --SelectComponents 0"
+UNPACKER+=" --ReqMode $_reqmod"
+UNPACKER+=" --ReqTint $_reqtint"
+UNPACKER+=" --PulserMode $_pulmode"
+UNPACKER+=" --PulMulMin $_pulmulmin"
+UNPACKER+=" --PulTotMin $_pultotmin"
+UNPACKER+=" --PulTotMax $_pultotmax"
+UNPACKER+=" --ToffTof $_tofftof"
+UNPACKER+=" --RefModType -1"
+UNPACKER+=" --RefModId    0"
+UNPACKER+=" --RefCtrType  0"
+UNPACKER+=" --RefCtrId    0"
+if  [[ $_batch = 1 ]]; then 
+UNPACKER+=" --control static"
+UNPACKER+=" --log-to-file Unp$_iUnp.out"
+nohup /home/cbm/software/cbmroot/git/cbmroot/build/bin/MQ/unpacker/$UNPACKER &
+else 
+xterm -geometry 110x23+520+0 -hold -e /home/cbm/software/cbmroot/git/cbmroot/build/bin/MQ/unpacker/$UNPACKER & 
+fi
+
+(( _iUnp -= 1 ))
+done 
+
+HITBUILDER="HitBuilderTof"
+HITBUILDER+=" --id hitbuilder1"
+HITBUILDER+=" --channel-config name=tofdigis,type=pull,method=bind,rateLogging=0,transport=zeromq,address=tcp://*:5656"
+HITBUILDER+=" --channel-config name=parameters,type=req,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:5005"
+HITBUILDER+=" --channel-config name=syscmd,type=sub,method=connect,rateLogging=0,transport=zeromq,address=tcp://127.0.0.1:5666"
+#HITBUILDER+=" --channel-config name=tofhits,type=push,method=bind,transport=shmem,address=tcp://127.0.0.1:5557"
+#HITBUILDER+=" --channel-config name=tofcalib,type=push,method=bind,transport=shmem,address=tcp://127.0.0.1:5558"
+#HITBUILDER+=" --transport shmem"
+#HITBUILDER+=" --severity DEBUG"
+HITBUILDER+=" --severity INFO"
+#HITBUILDER+=" --severity WARN"
+HITBUILDER+=" --OutRootFile $_outdir$_digifile"
+#HITBUILDER+=" --MaxEvent 10000000"
+#HITBUILDER+=" --RunId 1581312162"   # v20a
+HITBUILDER+=" --RunId 1597162455"   # v20b
+HITBUILDER+=" --PulserMode $_pulmode"
+HITBUILDER+=" --PulMulMin $_pulmulmin"
+HITBUILDER+=" --PulTotMin $_pultotmin"
+HITBUILDER+=" --PulTotMax $_pultotmax"
+HITBUILDER+=" --PulDetRef $_puldetref"
+HITBUILDER+=" --ReqTint   $_reqtint"
+#HITBUILDER+=" --ReqBeam      20486" # diamond -> 0x00005006
+HITBUILDER+=" --BRefType  9"
+HITBUILDER+=" --BRefSm    1"
+HITBUILDER+=" --BRefDet   0"
+HITBUILDER+=" --DutType   9"
+HITBUILDER+=" --DutSm     1"
+HITBUILDER+=" --DutRpc    0"
+HITBUILDER+=" --SelType   9"
+HITBUILDER+=" --SelSm     0"
+HITBUILDER+=" --SelRpc    1"
+HITBUILDER+=" --Sel2Type  9"
+HITBUILDER+=" --Sel2Sm    1"
+HITBUILDER+=" --Sel2Rpc   1"
+
+if  [[ $_batch = 1 ]]; then 
+HITBUILDER+=" --control static"
+HITBUILDER+=" --log-to-file HitBuild.out"
+nohup  /home/cbm/software/cbmroot/git/cbmroot/build/bin/MQ/hitbuilder/$HITBUILDER &
+else
+xterm -geometry 120x23+1400+0 -hold -e /home/cbm/software/cbmroot/git/cbmroot/build/bin/MQ/hitbuilder/$HITBUILDER &
+fi
diff --git a/macro/beamtime/fit_yPos.C b/macro/beamtime/fit_yPos.C
new file mode 100644
index 0000000000000000000000000000000000000000..0b8205fabdf1d579e8dc759a842f4d59921646b3
--- /dev/null
+++ b/macro/beamtime/fit_yPos.C
@@ -0,0 +1,56 @@
+void fit_yPos(Int_t SmT = 0, Int_t iSm = 0, Int_t iRpc = 0) {
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2);
+  TCanvas* can = new TCanvas("can", "can", 50, 0, 800, 800);
+  can->Divide(1, 2);
+
+  gPad->SetFillColor(0);
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(kTRUE);
+  gStyle->SetOptFit(111);
+
+  gROOT->cd();
+  gROOT->SetDirLevel(1);
+
+  TH1* h;
+  TH1* h1;
+  TH2* h2;
+  TH1* h2y;
+  void fit_ybox(const char* hname);
+  void fit_ybox(const char* hname, Double_t ysize);
+  // if (h!=NULL) h->Delete();
+
+  can->cd(1);
+  gROOT->cd();
+  gROOT->LoadMacro("fit_ybox.h");
+   ROOT::Math::Minimizer* minimum =
+      ROOT::Math::Factory::CreateMinimizer("Minuit", "Migrad");
+  minimum->SetMaxFunctionCalls(100000);
+  minimum->SetTolerance(0.1);
+  minimum->SetPrintLevel(3);
+
+  TFitter::SetMaxIterations(500000);
+  TFitter::SetPrecision(0.1);
+  TFitter::SetDefaultFitter("Minuit2");
+
+  TString hname2 = Form("cl_SmT%d_sm%03d_rpc%03d_Pos", SmT, iSm, iRpc);
+  h2             = (TH2*) gROOT->FindObjectAny(hname2);
+  if (h2 != NULL) {
+    h2->Draw("colz");
+    gPad->SetLogz();
+    h2->ProfileX()->Draw("same");
+
+    can->cd(2);
+    h2y = h2->ProjectionY();
+    cout << " Fit with ybox " << h2y->GetName() << endl;
+    fit_ybox((const char*) (h2y->GetName()));
+    if( 0) { //NULL != gMinuit ) {
+      cout << "Minuit ended with " << gMinuit->fCstatu<<endl;
+    }
+  } else {
+    cout << hname2 << " not found" << endl;
+  }
+
+
+  can->SaveAs(Form("Ypos%01d_%01d_%01d.pdf", SmT, iSm, iRpc));
+}
diff --git a/macro/beamtime/mcbm2018/ana_digi_cal_all.C b/macro/beamtime/mcbm2018/ana_digi_cal_all.C
new file mode 100644
index 0000000000000000000000000000000000000000..83c27165a0218ee8faba8c1385b8097ac3f17cbc
--- /dev/null
+++ b/macro/beamtime/mcbm2018/ana_digi_cal_all.C
@@ -0,0 +1,454 @@
+void ana_digi_cal_all(Int_t nEvents      = 10000000,
+                      Int_t calMode      = 53,
+                      Int_t calSel       = 0,
+                      Int_t calSm        = 900,
+                      Int_t RefSel       = 1,
+                      TString cFileId    = "Test",
+                      Int_t iCalSet      = 910601600,
+                      Bool_t bOut        = 0,
+                      Int_t iSel2        = 0,
+                      Double_t dDeadtime = 50,
+                      TString cCalId     = "XXX",
+                      Int_t iPlot        = 1) {
+  Int_t iVerbose = 1;
+  Int_t iBugCor  = 0;
+  //Specify log level (INFO, DEBUG, DEBUG1, ...)
+  //TString logLevel = "FATAL";
+  //TString logLevel = "ERROR";
+  TString logLevel = "INFO";
+  //TString logLevel = "DEBUG";
+  //TString logLevel = "DEBUG1";
+  //TString logLevel = "DEBUG2";
+  //TString logLevel = "DEBUG3";
+  FairLogger::GetLogger();
+  gLogger->SetLogScreenLevel(logLevel);
+  gLogger->SetLogVerbosityLevel("VERYHIGH");
+
+  TString workDir = gSystem->Getenv("VMCWORKDIR");
+  /*
+   TString workDir    = (TString)gInterpreter->ProcessLine(".! pwd");
+   cout << "workdir = "<< workDir.Data() << endl;
+   return;
+  */
+  TString paramDir  = workDir + "/macro/beamtime/mcbm2018/";
+  TString ParFile   = paramDir + "data/" + cFileId + ".params.root";
+  TString InputFile = paramDir + "data/" + cFileId + ".root";
+  // TString InputFile  =  "./data/" + cFileId + ".root";
+  TString OutputFile =
+    paramDir + "data/digidev_" + cFileId
+    + Form("_%09d_%03d_%02.0f_Cal", iCalSet, iSel2, dDeadtime) + cCalId
+    + ".out.root";
+
+  TString shcmd = "rm -v " + ParFile;
+  gSystem->Exec(shcmd.Data());
+
+  TList* parFileList = new TList();
+
+  TString FId    = cFileId;
+  TString TofGeo = "v18m_mcbm";
+  cout << "Geometry version " << TofGeo << endl;
+/*
+  TObjString* tofDigiFile = new TObjString(
+    workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par");  // TOF digi file
+  parFileList->Add(tofDigiFile);
+*/
+  //   TObjString tofDigiBdfFile = new TObjString( paramDir + "/tof." + FPar + "digibdf.par");
+  TObjString* tofDigiBdfFile =
+    new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digibdf.par");
+  parFileList->Add(tofDigiBdfFile);
+
+  TString geoDir  = gSystem->Getenv("VMCWORKDIR");
+  TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+  TFile* fgeo     = new TFile(geoFile);
+  TGeoManager* geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+  if (NULL == geoMan) {
+    cout << "<E> FAIRGeom not found in geoFile" << endl;
+    return;
+  }
+
+  if (0) {
+    TGeoVolume* master = geoMan->GetTopVolume();
+    master->SetVisContainers(1);
+    master->Draw("ogl");
+  }
+
+  // -----   Reconstruction run   -------------------------------------------
+  FairRunAna* run = new FairRunAna();
+  run->SetInputFile(InputFile.Data());
+  //run->AddFriend(InputFile.Data());
+  // run->SetOutputFile(OutputFile);
+  //run->SetSink( new FairRootFileSink( OutputFile.Data() ) );
+  run->SetUserOutputFileName(OutputFile.Data());
+  run->SetSink(new FairRootFileSink(run->GetUserOutputFileName()));
+  CbmTofEventClusterizer* tofClust =
+    new CbmTofEventClusterizer("TOF Event Clusterizer", iVerbose, bOut);
+
+  tofClust->SetCalMode(calMode);
+  tofClust->SetCalSel(calSel);
+  tofClust->SetCaldXdYMax(3.);  // geometrical matching window in cm
+  tofClust->SetCalCluMulMax(
+    5.);  // Max Counter Cluster Multiplicity for filling calib histos
+  tofClust->SetCalRpc(calSm);   // select detector for calibration update
+  tofClust->SetTRefId(RefSel);  // reference trigger for offset calculation
+  tofClust->SetTotMax(20.);     // Tot upper limit for walk corection
+  tofClust->SetTotMin(
+    0.01);  //(12000.);  // Tot lower limit for walk correction
+  tofClust->SetTotPreRange(
+    5.);  // effective lower Tot limit  in ns from peak position
+  tofClust->SetTotMean(5.);       // Tot calibration target value in ns
+  tofClust->SetMaxTimeDist(1.0);  // default cluster range in ns
+  //tofClust->SetMaxTimeDist(0.);       //Deb// default cluster range in ns
+  tofClust->SetDelTofMax(
+    5.);  // acceptance range for cluster distance in ns (!)
+  tofClust->SetSel2MulMax(3);  // limit Multiplicity in 2nd selector
+  tofClust->SetChannelDeadtime(dDeadtime);  // artificial deadtime in ns
+  tofClust->SetEnableAvWalk(kFALSE);
+  //tofClust->SetEnableMatchPosScaling(kFALSE); // turn off projection to nominal target
+  tofClust->SetYFitMin(1.E4);
+  tofClust->SetToDAv(0.04);
+  // tofClust->SetTimePeriod(25600.);       // ignore coarse time
+  // tofClust->SetCorMode(iBugCor);         // correct missing hits
+  //tofClust->SetIdMode(0);                  // calibrate on counter level
+  tofClust->SetIdMode(1);  // calibrate on module level
+  //   tofClust->SetDeadStrips(15,23);   // declare dead strip for T0M3,Rpc0,Strip 23
+  //tofClust->SetDeadStrips(25,16);   // declare non-existant diamond strip (#5) dead
+
+  Int_t calSelRead = calSel;
+  if (calSel < 0) calSelRead = 0;
+  TString cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                        cFileId.Data(),
+                        iCalSet,
+                        calMode,
+                        calSelRead);
+  if (cCalId != "XXX")
+    cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                  cCalId.Data(),
+                  iCalSet,
+                  calMode,
+                  calSelRead);
+  tofClust->SetCalParFileName(cFname);
+  TString cOutFname =
+    Form("tofClust_%s_set%09d.hst.root", cFileId.Data(), iCalSet);
+  tofClust->SetOutHstFileName(cOutFname);
+
+  TString cAnaFile =
+    Form("%s_%09d%03d_tofAna.hst.root", cFileId.Data(), iCalSet, iSel2);
+
+  switch (calMode) {
+    case -1:                      // initial check of raw data
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      //tofClust->SetTotMin(1.);
+      tofClust->SetTRefDifMax(26000.);  // in ns
+      tofClust->PosYMaxScal(10000.);    // in % of length
+      tofClust->SetMaxTimeDist(0.);     // no cluster building
+      //tofClust->SetTimePeriod(25600.);       // inspect coarse time
+      break;
+    case 0:                       // initial calibration
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      //tofClust->SetTotMin(1.);
+      tofClust->SetTRefDifMax(1000.);  // in ns
+      tofClust->PosYMaxScal(10.);      // in % of length
+      tofClust->SetMaxTimeDist(0.);    // no cluster building
+      break;
+    case 1:                       // save offsets, update walks, for diamonds
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      tofClust->SetTRefDifMax(6.25);  // in ns
+      //tofClust->SetTimePeriod(6.25);       // inspect coarse time
+      tofClust->PosYMaxScal(10.);  // in % of length
+      break;
+    case 11:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(3.0);   // in % of length
+      break;
+    case 21:
+      tofClust->SetTRefDifMax(3.0);  // in ns
+      tofClust->PosYMaxScal(2.0);    // in % of length
+      break;
+    case 31:
+      tofClust->SetTRefDifMax(3.);  // in ns
+      tofClust->PosYMaxScal(1.);    // in % of length
+      break;
+    case 41:
+      tofClust->SetTRefDifMax(2.0);  // in ns
+      tofClust->PosYMaxScal(0.9);    // in % of length
+      break;
+    case 51:
+      tofClust->SetTRefDifMax(2.0);  // in ns
+      tofClust->PosYMaxScal(0.8);    // in % of length
+      break;
+    case 61:
+      tofClust->SetTRefDifMax(1.5);  // in ns
+      tofClust->PosYMaxScal(0.75);   // in % of length
+      break;
+    case 71:
+      tofClust->SetTRefDifMax(0.8);  // in ns
+      tofClust->PosYMaxScal(0.6);    // in % of length
+      break;
+
+    case 2:                           // time difference calibration
+      tofClust->SetTRefDifMax(300.);  // in ns
+      tofClust->PosYMaxScal(1000.);   //in % of length
+      break;
+
+    case 3:                           // time offsets
+      tofClust->SetTRefDifMax(200.);  // in ns
+      tofClust->PosYMaxScal(100.);    //in % of length
+      tofClust->SetMaxTimeDist(0.);   // no cluster building
+      break;
+    case 12:
+    case 13:
+      tofClust->SetTRefDifMax(100.);  // in ns
+      tofClust->PosYMaxScal(10.);     //in % of length
+      break;
+    case 22:
+    case 23:
+      tofClust->SetTRefDifMax(50.);  // in ns
+      tofClust->PosYMaxScal(5.);     //in % of length
+      break;
+    case 32:
+    case 33:
+      tofClust->SetTRefDifMax(25.);  // in ns
+      tofClust->PosYMaxScal(4.);     //in % of length
+      break;
+    case 42:
+    case 43:
+      tofClust->SetTRefDifMax(12.);  // in ns
+      tofClust->PosYMaxScal(2.);     //in % of length
+      break;
+    case 52:
+    case 53:
+      tofClust->SetTRefDifMax(10.);  // in ns
+      tofClust->PosYMaxScal(1.5);    //in % of length
+      break;
+    case 62:
+    case 63:
+      tofClust->SetTRefDifMax(10.);  // in ns
+      tofClust->PosYMaxScal(1.);     //in % of length
+      break;
+    case 72:
+    case 73:
+      tofClust->SetTRefDifMax(10.);  // in ns
+      tofClust->PosYMaxScal(0.9);    //in % of length
+      break;
+    case 82:
+    case 83:
+      tofClust->SetTRefDifMax(10.);  // in ns
+      tofClust->PosYMaxScal(0.8);    //in % of length
+      break;
+    case 92:
+    case 93:
+      tofClust->SetTRefDifMax(10.);  // in ns
+      tofClust->PosYMaxScal(0.75);   //in % of length
+      break;
+
+    case 4:  // velocity dependence (DelTOF)
+    case 14:
+      tofClust->SetTRefDifMax(25.);  // in ns
+      tofClust->PosYMaxScal(2.0);    //in % of length
+      break;
+    case 24:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(1.5);   //in % of length
+      break;
+    case 34:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(1.2);   //in % of length
+      break;
+    case 44:
+      tofClust->SetTRefDifMax(3.5);  // in ns
+      tofClust->PosYMaxScal(1.0);    //in % of length
+      break;
+    case 54:
+      tofClust->SetTRefDifMax(3.0);  // in ns
+      tofClust->PosYMaxScal(0.9);    //in % of length
+      break;
+    case 64:
+      tofClust->SetTRefDifMax(2.5);  // in ns
+      tofClust->PosYMaxScal(0.8);    //in % of length
+      break;
+    case 74:
+      tofClust->SetTRefDifMax(2.0);  // in ns
+      tofClust->PosYMaxScal(0.7);    //in % of length
+      break;
+    default:
+      cout << "<E> Calib mode not implemented! stop execution of script"
+           << endl;
+      return;
+  }
+
+  Int_t iBRef    = iCalSet % 1000;
+  Int_t iSet     = (iCalSet - iBRef) / 1000;
+  Int_t iRSel    = 0;
+  Int_t iRSelTyp = 0;
+  Int_t iRSelSm  = 0;
+  Int_t iRSelRpc = 0;
+  iRSel          = iBRef;  // use diamond
+
+  Int_t iRSelin = iRSel;
+  iRSelRpc      = iRSel % 10;
+  iRSelTyp      = (iRSel - iRSelRpc) / 10;
+  iRSelSm       = iRSelTyp % 10;
+  iRSelTyp      = (iRSelTyp - iRSelSm) / 10;
+
+  tofClust->SetBeamRefId(iRSelTyp);  // define Beam reference counter
+  tofClust->SetBeamRefSm(iRSelSm);
+  tofClust->SetBeamRefDet(iRSelRpc);
+  tofClust->SetBeamAddRefMul(-1);
+  tofClust->SetBeamRefMulMax(3);
+
+  Int_t iSel2in  = iSel2;
+  Int_t iSel2Rpc = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Rpc) / 10;
+  Int_t iSel2Sm  = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Sm) / 10;
+
+  tofClust->SetSel2Id(iSel2);
+  tofClust->SetSel2Sm(iSel2Sm);
+  tofClust->SetSel2Rpc(iSel2Rpc);
+
+  Int_t iRef    = iSet % 1000;
+  Int_t iDut    = (iSet - iRef) / 1000;
+  Int_t iDutRpc = iDut % 10;
+  iDut          = (iDut - iDutRpc) / 10;
+  Int_t iDutSm  = iDut % 10;
+  iDut          = (iDut - iDutSm) / 10;
+
+  tofClust->SetDutId(iDut);
+  tofClust->SetDutSm(iDutSm);
+  tofClust->SetDutRpc(iDutRpc);
+
+  Int_t iRefRpc = iRef % 10;
+  iRef          = (iRef - iRefRpc) / 10;
+  Int_t iRefSm  = iRef % 10;
+  iRef          = (iRef - iRefSm) / 10;
+
+  tofClust->SetSelId(iRef);
+  tofClust->SetSelSm(iRefSm);
+  tofClust->SetSelRpc(iRefRpc);
+
+  run->AddTask(tofClust);
+
+  cout << "Run with iRSel = " << iRSel << ", iSel2 = " << iSel2in << endl;
+
+
+  // -----  Parameter database   --------------------------------------------
+
+  FairRuntimeDb* rtdb       = run->GetRuntimeDb();
+  Bool_t kParameterMerged   = kTRUE;
+  FairParRootFileIo* parIo2 = new FairParRootFileIo(kParameterMerged);
+  parIo2->open(ParFile.Data(), "UPDATE");
+  parIo2->print();
+  rtdb->setFirstInput(parIo2);
+
+  FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
+  parIo1->open(parFileList, "in");
+  parIo1->print();
+  rtdb->setSecondInput(parIo1);
+  rtdb->print();
+  rtdb->printParamContexts();
+
+  //  FairParRootFileIo* parInput1 = new FairParRootFileIo();
+  //  parInput1->open(ParFile.Data());
+  //  rtdb->setFirstInput(parInput1);
+
+  // -----   Intialise and run   --------------------------------------------
+  run->Init();
+  cout << "Starting run" << endl;
+  run->Run(0, nEvents);
+  //tofClust->Finish();
+  // ------------------------------------------------------------------------
+  // default display
+  /*
+  TString Display_Status = "pl_over_Mat04D4best.C";
+  TString Display_Funct = "pl_over_Mat04D4best()";  
+  gROOT->LoadMacro(Display_Status);
+  */
+
+  gROOT->LoadMacro("save_hst.C");
+
+  gROOT->LoadMacro("fit_ybox.h");
+  gROOT->LoadMacro("pl_all_CluMul.C");
+  gROOT->LoadMacro("pl_all_CluRate.C");
+  gROOT->LoadMacro("pl_all_CluPosEvol.C");
+  gROOT->LoadMacro("pl_all_CluTimeEvol.C");
+  gROOT->LoadMacro("pl_over_cluSel.C");
+  gROOT->LoadMacro("pl_over_clu.C");
+  gROOT->LoadMacro("pl_over_Walk2.C");
+  gROOT->LoadMacro("pl_all_dTSel.C");
+  gROOT->LoadMacro("pl_over_MatD4sel.C");
+  gROOT->LoadMacro("pl_all_Sel2D.C");
+  gROOT->LoadMacro("pl_all_2D.C");
+
+  if (iPlot) {
+
+    switch (iSet) {
+      default:
+        for (Int_t iOpt = 0; iOpt < 8; iOpt++) {
+          for (Int_t iSel = 0; iSel < 2; iSel++) {
+            gInterpreter->ProcessLine(Form("pl_all_Sel2D(%d,%d)", iOpt, iSel));
+          }
+        }
+
+        for (Int_t iOpt = 0; iOpt < 12; iOpt++) {
+          gInterpreter->ProcessLine(Form("pl_all_2D(%d)", iOpt));
+        }
+        /*
+	gInterpreter->ProcessLine("pl_over_clu(0,0,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,4)");
+	
+	gInterpreter->ProcessLine("pl_over_clu(5,0,0)");
+	gInterpreter->ProcessLine("pl_over_cluSel(0,5,0,0)");
+	gInterpreter->ProcessLine("pl_over_cluSel(1,5,0,0)");
+	
+	for(Int_t iSm=0; iSm<3; iSm++)
+	for (Int_t iRpc=0; iRpc<5; iRpc++)
+	for (Int_t iSel=0; iSel<2; iSel++){
+	gInterpreter->ProcessLine(Form("pl_over_cluSel(%d,0,%d,%d)",iSel,iSm,iRpc));
+	gInterpreter->ProcessLine(Form("pl_over_Walk2(%d,0,%d,%d)",iSel,iSm,iRpc));
+	}
+      */
+        gInterpreter->ProcessLine("pl_all_CluMul()");
+        gInterpreter->ProcessLine("pl_all_CluRate()");
+        gInterpreter->ProcessLine("pl_all_CluRate(5,1)");
+        gInterpreter->ProcessLine("pl_all_CluPosEvol()");
+        gInterpreter->ProcessLine("pl_all_CluTimeEvol()");
+        gInterpreter->ProcessLine("pl_all_dTSel()");
+
+        //  gInterpreter->ProcessLine("pl_over_MatD4sel()");
+        //  gInterpreter->ProcessLine(Display_Funct.Data());
+        break;
+        ;
+    }
+  }
+  TString FSave = Form("save_hst(\"CluStatus%d_%d_Cal_%s.hst.root\")",
+                       iCalSet,
+                       iSel2in,
+                       cCalId.Data());
+  gInterpreter->ProcessLine(FSave.Data());
+}
diff --git a/macro/beamtime/mcbm2018/ana_digi_cali.C b/macro/beamtime/mcbm2018/ana_digi_cali.C
new file mode 100644
index 0000000000000000000000000000000000000000..6fa1fa55ae60e48ce64f43a0b0fa5153b1106969
--- /dev/null
+++ b/macro/beamtime/mcbm2018/ana_digi_cali.C
@@ -0,0 +1,460 @@
+void ana_digi_cali(Int_t nEvents      = 10000000,
+                  Int_t calMode      = 53,
+                  Int_t calSel       = 0,
+                  Int_t calSm        = 900,
+                  Int_t RefSel       = 1,
+                  TString cFileId    = "Test",
+                  Int_t iCalSet      = 910601600,
+                  Bool_t bOut        = 0,
+                  Int_t iSel2        = 0,
+                  Double_t dDeadtime = 50,
+                  TString cCalId     = "XXX",
+                  Int_t iPlot        = 1) {
+  Int_t iVerbose = 1;
+  Int_t iBugCor  = 0;
+  //Specify log level (INFO, DEBUG, DEBUG1, ...)
+  //TString logLevel = "FATAL";
+  //TString logLevel = "ERROR";
+  TString logLevel = "INFO";
+  //TString logLevel = "DEBUG";
+  //TString logLevel = "DEBUG1";
+  //TString logLevel = "DEBUG2";
+  //TString logLevel = "DEBUG3";
+  FairLogger::GetLogger();
+  gLogger->SetLogScreenLevel(logLevel);
+  gLogger->SetLogVerbosityLevel("MEDIUM");
+
+  TString workDir = gSystem->Getenv("VMCWORKDIR");
+  /*
+   TString workDir    = (TString)gInterpreter->ProcessLine(".! pwd");
+   cout << "workdir = "<< workDir.Data() << endl;
+   return;
+  */
+  //TString paramDir  = workDir + "/macro/beamtime/mcbm2018/";
+  TString paramDir  = "./"; 
+  TString ParFile   = paramDir + "data/" + cFileId + ".params.root";
+  TString InputFile = paramDir + "data/" + cFileId + ".root";
+  TString OutputFile =
+    paramDir + "data/digidev_" + cFileId
+    + Form("_%09d_%03d_%02.0f_Cal", iCalSet, iSel2, dDeadtime) + cCalId
+    + ".out.root";
+
+  TString shcmd = "rm -v " + ParFile;
+  gSystem->Exec(shcmd.Data());
+
+  TList* parFileList = new TList();
+
+  TString FId    = cFileId;
+  TString TofGeo = "v18m_mcbm";
+  cout << "Geometry version " << TofGeo << endl;
+/*
+  TObjString* tofDigiFile = new TObjString(
+    workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par");  // TOF digi file
+  parFileList->Add(tofDigiFile);
+*/
+  //   TObjString tofDigiBdfFile = new TObjString( paramDir + "/tof." + FPar + "digibdf.par");
+  TObjString* tofDigiBdfFile =
+    new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digibdf.par");
+  parFileList->Add(tofDigiBdfFile);
+
+  TString geoDir  = gSystem->Getenv("VMCWORKDIR");
+  TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+  TFile* fgeo     = new TFile(geoFile);
+  TGeoManager* geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+  if (NULL == geoMan) {
+    cout << "<E> FAIRGeom not found in geoFile" << endl;
+    return;
+  }
+
+  if (0) {
+    TGeoVolume* master = geoMan->GetTopVolume();
+    master->SetVisContainers(1);
+    master->Draw("ogl");
+  }
+
+  // -----   Reconstruction run   -------------------------------------------
+  FairRunAna* run = new FairRunAna();
+  run->SetInputFile(InputFile.Data());
+  //run->AddFriend(InputFile.Data());
+  // run->SetOutputFile(OutputFile);
+  //run->SetSink( new FairRootFileSink( OutputFile.Data() ) );
+  run->SetUserOutputFileName(OutputFile.Data());
+  run->SetSink(new FairRootFileSink(run->GetUserOutputFileName()));
+  CbmTofEventClusterizer* tofClust =
+    new CbmTofEventClusterizer("TOF Event Clusterizer", iVerbose, bOut);
+
+  tofClust->SetCalMode(calMode);
+  tofClust->SetCalSel(calSel);
+  tofClust->SetCaldXdYMax(10.);  // geometrical matching window in cm
+  tofClust->SetCalCluMulMax(
+    3.);  // Max Counter Cluster Multiplicity for filling calib histos
+  tofClust->SetCalRpc(calSm);   // select detector for calibration update
+  tofClust->SetTRefId(RefSel);  // reference trigger for offset calculation
+  tofClust->SetTotMax(20.);     // Tot upper limit for walk corection
+  tofClust->SetTotMin(
+    0.01);  //(12000.);  // Tot lower limit for walk correction
+  tofClust->SetTotPreRange(
+    5.);  // effective lower Tot limit  in ns from peak position
+  tofClust->SetTotMean(5.);       // Tot calibration target value in ns
+  tofClust->SetMaxTimeDist(1.0);  // default cluster range in ns
+  //tofClust->SetMaxTimeDist(0.);       //Deb// default cluster range in ns
+  tofClust->SetDelTofMax(
+    10.);  // acceptance range for cluster distance in ns (!)
+  tofClust->SetSel2MulMax(3);  // limit Multiplicity in 2nd selector
+  tofClust->SetChannelDeadtime(dDeadtime);  // artificial deadtime in ns
+  tofClust->SetEnableAvWalk(kFALSE);
+  //tofClust->SetEnableMatchPosScaling(kFALSE); // turn off projection to nominal target
+  tofClust->SetYFitMin(1.E4);
+  // tofClust->SetTimePeriod(25600.);       // ignore coarse time
+  // tofClust->SetCorMode(iBugCor);         // correct missing hits
+  tofClust->SetIdMode(1);  // calibrate on module level
+  //   tofClust->SetDeadStrips(15,23);   // declare dead strip for T0M3,Rpc0,Strip 23
+  //   tofClust->SetDeadStrips(27,30);   // declare dead strip for T6M0,Rpc1,Strip 30
+  //   tofClust->SetDeadStrips(28,19);   // declare dead strip for T7M0,Rpc0,Strip 19
+  //tofClust->SetDeadStrips(25,16);   // declare non-existant diamond strip (#5) dead
+
+  Int_t calSelRead = calSel;
+  if (calSel < 0) calSelRead = 0;
+  TString cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                        cFileId.Data(),
+                        iCalSet,
+                        calMode,
+                        calSelRead);
+  if (cCalId != "XXX")
+    cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                  cCalId.Data(),
+                  iCalSet,
+                  calMode,
+                  calSelRead);
+  tofClust->SetCalParFileName(cFname);
+  TString cOutFname =
+    Form("tofClust_%s_set%09d.hst.root", cFileId.Data(), iCalSet);
+  tofClust->SetOutHstFileName(cOutFname);
+
+  TString cAnaFile =
+    Form("%s_%09d%03d_tofAna.hst.root", cFileId.Data(), iCalSet, iSel2);
+
+  switch (calMode) {
+    case -1:                      // initial check of raw data
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      //tofClust->SetTotMin(1.);
+      tofClust->SetTRefDifMax(26000.);  // in ns
+      tofClust->PosYMaxScal(10000.);    // in % of length
+      tofClust->SetMaxTimeDist(0.);     // no cluster building
+      //tofClust->SetTimePeriod(25600.);       // inspect coarse time
+      break;
+    case 0:                       // initial calibration
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      //tofClust->SetTotMin(1.);
+      tofClust->SetTRefDifMax(1000.);  // in ns
+      tofClust->PosYMaxScal(10.);      // in % of length
+      tofClust->SetMaxTimeDist(0.);    // no cluster building
+      break;
+    case 1:                       // save offsets, update walks, for diamonds
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      tofClust->SetTRefDifMax(6.25);  // in ns
+      //tofClust->SetTimePeriod(6.25);       // inspect coarse time
+      tofClust->PosYMaxScal(10.);  // in % of length
+      break;
+    case 11:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(3.0);   // in % of length
+      break;
+    case 21:
+      tofClust->SetTRefDifMax(2.5);  // in ns
+      tofClust->PosYMaxScal(2.0);    // in % of length
+      break;
+    case 31:
+      tofClust->SetTRefDifMax(2.);  // in ns
+      tofClust->PosYMaxScal(1.5);   // in % of length
+      break;
+    case 41:
+      tofClust->SetTRefDifMax(1.5);  // in ns
+      tofClust->PosYMaxScal(0.8);    // in % of length
+      break;
+    case 51:
+      tofClust->SetTRefDifMax(1.2);  // in ns
+      tofClust->PosYMaxScal(0.7);    // in % of length
+      break;
+    case 61:
+      tofClust->SetTRefDifMax(1.0);  // in ns
+      tofClust->PosYMaxScal(0.75);   // in % of length
+      break;
+    case 71:
+      tofClust->SetTRefDifMax(0.8);  // in ns
+      tofClust->PosYMaxScal(0.6);    // in % of length
+      break;
+
+    case 2:                           // time difference calibration
+      tofClust->SetTRefDifMax(300.);  // in ns
+      tofClust->PosYMaxScal(1000.);   //in % of length
+      break;
+
+    case 3:                           // time offsets
+      tofClust->SetTRefDifMax(200.);  // in ns
+      tofClust->PosYMaxScal(100.);    //in % of length
+      tofClust->SetMaxTimeDist(0.);   // no cluster building
+      break;
+    case 12:
+    case 13:
+      tofClust->SetTRefDifMax(100.);  // in ns
+      tofClust->PosYMaxScal(10.);     //in % of length
+      break;
+    case 22:
+    case 23:
+      tofClust->SetTRefDifMax(50.);  // in ns
+      tofClust->PosYMaxScal(5.);     //in % of length
+      break;
+    case 32:
+    case 33:
+      tofClust->SetTRefDifMax(25.);  // in ns
+      tofClust->PosYMaxScal(4.);     //in % of length
+      break;
+    case 42:
+    case 43:
+      tofClust->SetTRefDifMax(12.);  // in ns
+      tofClust->PosYMaxScal(2.);     //in % of length
+      break;
+    case 52:
+    case 53:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(1.5);   //in % of length
+      break;
+    case 62:
+    case 63:
+      tofClust->SetTRefDifMax(3.);  // in ns
+      tofClust->PosYMaxScal(1.);    //in % of length
+      break;
+    case 72:
+    case 73:
+      tofClust->SetTRefDifMax(2.5);  // in ns
+      tofClust->PosYMaxScal(0.9);    //in % of length
+      break;
+    case 82:
+    case 83:
+      tofClust->SetTRefDifMax(2.0);  // in ns
+      tofClust->PosYMaxScal(0.8);    //in % of length
+      break;
+    case 92:
+    case 93:
+      tofClust->SetTRefDifMax(1.5);  // in ns
+      tofClust->PosYMaxScal(0.75);   //in % of length
+      break;
+
+    case 4:  // velocity dependence (DelTOF)
+    case 14:
+      tofClust->SetTRefDifMax(25.);  // in ns
+      tofClust->PosYMaxScal(2.0);    //in % of length
+      break;
+    case 24:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(1.5);   //in % of length
+      break;
+    case 34:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(1.2);   //in % of length
+      break;
+    case 44:
+      tofClust->SetTRefDifMax(3.);  // in ns
+      tofClust->PosYMaxScal(1.0);   //in % of length
+      break;
+    case 54:
+      tofClust->SetTRefDifMax(2.5);  // in ns
+      tofClust->PosYMaxScal(0.9);    //in % of length
+      break;
+    case 64:
+      tofClust->SetTRefDifMax(2.5);  // in ns
+      tofClust->PosYMaxScal(0.8);    //in % of length
+      break;
+    case 74:
+      tofClust->SetTRefDifMax(1.5);  // in ns
+      tofClust->PosYMaxScal(0.7);    //in % of length
+      break;
+    default:
+      cout << "<E> Calib mode not implemented! stop execution of script"
+           << endl;
+      return;
+  }
+
+  run->AddTask(tofClust);
+
+  Int_t iBRef    = iCalSet % 1000;
+  Int_t iSet     = (iCalSet - iBRef) / 1000;
+  Int_t iRSel    = 0;
+  Int_t iRSelTyp = 0;
+  Int_t iRSelSm  = 0;
+  Int_t iRSelRpc = 0;
+  iRSel          = iBRef;  // use diamond
+  if (iSel2 == 0) {
+    // iSel2=iBRef;
+  } else {
+    if (iSel2 < 0) iSel2 = -iSel2;
+  }
+
+  Int_t iRSelin = iRSel;
+  iRSelRpc      = iRSel % 10;
+  iRSelTyp      = (iRSel - iRSelRpc) / 10;
+  iRSelSm       = iRSelTyp % 10;
+  iRSelTyp      = (iRSelTyp - iRSelSm) / 10;
+
+  tofClust->SetBeamRefId(iRSelTyp);  // define Beam reference counter
+  tofClust->SetBeamRefSm(iRSelSm);
+  tofClust->SetBeamRefDet(iRSelRpc);
+  tofClust->SetBeamAddRefMul(-1);
+  tofClust->SetBeamRefMulMax(3);
+
+  Int_t iSel2in  = iSel2;
+  Int_t iSel2Rpc = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Rpc) / 10;
+  Int_t iSel2Sm  = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Sm) / 10;
+  if (iSel2 > 0) {
+    tofClust->SetSel2Id(iSel2);
+    tofClust->SetSel2Sm(iSel2Sm);
+    tofClust->SetSel2Rpc(iSel2Rpc);
+  }
+
+  Int_t iRef    = iSet % 1000;
+  Int_t iDut    = (iSet - iRef) / 1000;
+  Int_t iDutRpc = iDut % 10;
+  iDut          = (iDut - iDutRpc) / 10;
+  Int_t iDutSm  = iDut % 10;
+  iDut          = (iDut - iDutSm) / 10;
+
+  tofClust->SetDutId(iDut);
+  tofClust->SetDutSm(iDutSm);
+  tofClust->SetDutRpc(iDutRpc);
+
+  Int_t iRefRpc = iRef % 10;
+  iRef          = (iRef - iRefRpc) / 10;
+  Int_t iRefSm  = iRef % 10;
+  iRef          = (iRef - iRefSm) / 10;
+
+  tofClust->SetSelId(iRef);
+  tofClust->SetSelSm(iRefSm);
+  tofClust->SetSelRpc(iRefRpc);
+
+  cout << "Run with iRSel = " << iRSel << ", iSel2 = " << iSel2in << endl;
+
+
+  // -----  Parameter database   --------------------------------------------
+
+  FairRuntimeDb* rtdb       = run->GetRuntimeDb();
+  Bool_t kParameterMerged   = kTRUE;
+  FairParRootFileIo* parIo2 = new FairParRootFileIo(kParameterMerged);
+  parIo2->open(ParFile.Data(), "UPDATE");
+  parIo2->print();
+  rtdb->setFirstInput(parIo2);
+
+  FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
+  parIo1->open(parFileList, "in");
+  parIo1->print();
+  rtdb->setSecondInput(parIo1);
+  rtdb->print();
+  rtdb->printParamContexts();
+
+  //  FairParRootFileIo* parInput1 = new FairParRootFileIo();
+  //  parInput1->open(ParFile.Data());
+  //  rtdb->setFirstInput(parInput1);
+
+  // -----   Intialise and run   --------------------------------------------
+  run->Init();
+  cout << "Starting run" << endl;
+  run->Run(0, nEvents);
+  //tofClust->Finish();
+  // ------------------------------------------------------------------------
+  // default display
+  /*
+  TString Display_Status = "pl_over_Mat04D4best.C";
+  TString Display_Funct = "pl_over_Mat04D4best()";  
+  gROOT->LoadMacro(Display_Status);
+  */
+
+  gROOT->LoadMacro("save_hst.C");
+
+  gROOT->LoadMacro("fit_ybox.h");
+  gROOT->LoadMacro("pl_all_CluMul.C");
+  gROOT->LoadMacro("pl_all_CluRate.C");
+  gROOT->LoadMacro("pl_all_CluPosEvol.C");
+  gROOT->LoadMacro("pl_all_CluTimeEvol.C");
+  gROOT->LoadMacro("pl_over_cluSel.C");
+  gROOT->LoadMacro("pl_over_clu.C");
+  gROOT->LoadMacro("pl_over_Walk2.C");
+  gROOT->LoadMacro("pl_all_dTSel.C");
+  gROOT->LoadMacro("pl_over_MatD4sel.C");
+  gROOT->LoadMacro("pl_all_Sel2D.C");
+  gROOT->LoadMacro("pl_all_2D.C");
+
+  if (iPlot) {
+
+    switch (iSet) {
+      default:
+        for (Int_t iOpt = 0; iOpt < 7; iOpt++) {
+          for (Int_t iSel = 0; iSel < 2; iSel++) {
+            gInterpreter->ProcessLine(Form("pl_all_Sel2D(%d,%d)", iOpt, iSel));
+          }
+        }
+
+        for (Int_t iOpt = 0; iOpt < 12; iOpt++) {
+          gInterpreter->ProcessLine(Form("pl_all_2D(%d)", iOpt));
+        }
+        /*
+	gInterpreter->ProcessLine("pl_over_clu(0,0,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,0,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,1,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,2,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,3,4)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,0)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,1)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,2)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,3)");
+	gInterpreter->ProcessLine("pl_over_clu(0,4,4)");
+	
+	gInterpreter->ProcessLine("pl_over_clu(5,0,0)");
+	gInterpreter->ProcessLine("pl_over_cluSel(0,5,0,0)");
+	gInterpreter->ProcessLine("pl_over_cluSel(1,5,0,0)");
+	
+	for(Int_t iSm=0; iSm<3; iSm++)
+	for (Int_t iRpc=0; iRpc<5; iRpc++)
+	for (Int_t iSel=0; iSel<2; iSel++){
+	gInterpreter->ProcessLine(Form("pl_over_cluSel(%d,0,%d,%d)",iSel,iSm,iRpc));
+	gInterpreter->ProcessLine(Form("pl_over_Walk2(%d,0,%d,%d)",iSel,iSm,iRpc));
+	}
+      */
+        gInterpreter->ProcessLine("pl_all_CluMul()");
+        gInterpreter->ProcessLine("pl_all_CluRate()");
+        gInterpreter->ProcessLine("pl_all_CluRate(5,1)");
+        gInterpreter->ProcessLine("pl_all_CluPosEvol()");
+        gInterpreter->ProcessLine("pl_all_CluTimeEvol()");
+        gInterpreter->ProcessLine("pl_all_dTSel()");
+
+        //  gInterpreter->ProcessLine("pl_over_MatD4sel()");
+        //  gInterpreter->ProcessLine(Display_Funct.Data());
+        break;
+        ;
+    }
+  }
+  TString FSave = Form("save_hst(\"CluStatus%d_%d_Cal_%s.hst.root\")",
+                       iCalSet,
+                       iSel2in,
+                       cCalId.Data());
+  gInterpreter->ProcessLine(FSave.Data());
+}
diff --git a/macro/beamtime/mcbm2018/ana_trksi.C b/macro/beamtime/mcbm2018/ana_trksi.C
new file mode 100644
index 0000000000000000000000000000000000000000..2f2c9ecb2c96ec28eb7a84996207d7e09dc32039
--- /dev/null
+++ b/macro/beamtime/mcbm2018/ana_trksi.C
@@ -0,0 +1,709 @@
+void ana_trksi(Int_t nEvents        = 10000,
+              Int_t iSel           = 1,
+              Int_t iGenCor        = 1,
+              TString cFileId      = "48.50.7.1",
+              TString cSet         = "000010020",
+              Int_t iSel2          = 20,
+              Int_t iTrackingSetup = 2,
+              Double_t dScalFac    = 1.,
+              Double_t dChi2Lim2   = 500.,
+              Double_t dDeadtime   = 50,
+              TString cCalId       = "",
+              Int_t iAnaCor        = 1,
+              Bool_t bUseSigCalib  = kFALSE,
+              Int_t iCalSet        = 30040500,
+              Int_t iCalOpt        = 1,
+              Int_t iMc            = 0) {
+  Int_t iVerbose = 1;
+  if (cCalId == "") cCalId = cFileId;
+  TString FId = cFileId;
+  TString cRun(FId(0, 3));
+  Int_t iRun = cRun.Atoi();
+  // Specify log level (INFO, DEBUG, DEBUG1, ...)
+  //TString logLevel = "FATAL";
+  //TString logLevel = "ERROR";
+  TString logLevel = "INFO";
+  //TString logLevel = "DEBUG";
+  //TString logLevel = "DEBUG1";
+  //TString logLevel = "DEBUG2";
+  //TString logLevel = "DEBUG3";
+  TString workDir  = gSystem->Getenv("VMCWORKDIR");
+  //TString paramDir = workDir + "/macro/beamtime/mcbm2018";
+  TString paramDir       = ".";
+
+  TString ParFile       = paramDir + "/data/" + cFileId.Data() + ".params.root";
+  TString InputFile     = paramDir + "/data/" + cFileId.Data() + ".root";
+  TString InputDigiFile = paramDir + "/data/digidev_" + cFileId.Data()
+                          + Form("_%s_%02.0f_Cal", cSet.Data(), dDeadtime)
+                          + cCalId + ".out.root";
+  if (iMc == 1) {
+    InputFile     = paramDir + "/data/" + cFileId.Data() + ".raw.root";
+    InputDigiFile = paramDir + "/data/" + cFileId.Data() + ".rec.root";
+    iRun          = 700;
+  }
+  TString OutputFile = paramDir + "/data/hits_" + cFileId.Data()
+                       + Form("_%s_%06d_%03d", cSet.Data(), iSel, iSel2)
+                       + ".out.root";
+  TString cHstFile =
+    paramDir
+    + Form(
+      "/hst/%s_%03.0f_%s_%06d_%03d_%03.1f_%03.1f_trk%03d_Cal%s_Ana.hst.root",
+      cFileId.Data(),
+      dDeadtime,
+      cSet.Data(),
+      iSel,
+      iSel2,
+      dScalFac,
+      dChi2Lim2,
+      iTrackingSetup,
+      cCalId.Data());
+  TString cTrkFile = Form("%s_tofFindTracks.hst.root", cCalId.Data());
+  TString cAnaFile = Form("%s_TrkAnaTestBeam.hst.root", cFileId.Data());
+
+  cout << " InputDigiFile = " << InputDigiFile << endl;
+
+  TString shcmd = "rm -v " + ParFile;
+  gSystem->Exec(shcmd.Data());
+
+  TList* parFileList = new TList();
+  Int_t iGeo = 0;  //iMc;
+  if (iGeo == 0) {
+    TString TofGeo = "v18m_mcbm";  //default
+
+    cout << "Geometry version " << TofGeo << endl;
+
+    TObjString* tofDigiBdfFile = new TObjString(workDir + "/parameters/tof/tof_"
+                                                + TofGeo + ".digibdf.par");
+    parFileList->Add(tofDigiBdfFile);
+
+    TString geoDir  = gSystem->Getenv("VMCWORKDIR");
+    TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+    TFile* fgeo     = new TFile(geoFile);
+    TGeoManager* geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+    if (NULL == geoMan) {
+      cout << "<E> FAIRGeom not found in geoFile" << endl;
+      return;
+    }
+  } else {
+    TString setupName = "mcbm_beam_2019_03";
+    // -----   Load the geometry setup   -------------------------------------
+    TString setupFile =
+      workDir + "/geometry/setup/setup_" + setupName.Data() + ".C";
+    TString setupFunct = "setup_";
+    setupFunct         = setupFunct + setupName + "()";
+    std::cout << "-I- Loading macro " << setupFile << std::endl;
+    gROOT->LoadMacro(setupFile);
+    gROOT->ProcessLine(setupFunct);
+    CbmSetup* setup = CbmSetup::Instance();
+  }
+
+  // -----   Reconstruction run   -------------------------------------------
+  FairRunAna* run = new FairRunAna();
+  cout << "InputFile:     " << InputFile.Data() << endl;
+  cout << "InputDigiFile: " << InputDigiFile.Data() << endl;
+
+  //run->SetInputFile(InputFile.Data());
+  //run->AddFriend(InputDigiFile.Data());
+  run->SetInputFile(InputDigiFile.Data());
+  //run->AddFriend(InputFile.Data());
+  //run->SetOutputFile(OutputFile);
+  run->SetUserOutputFileName(OutputFile.Data());
+  run->SetSink(new FairRootFileSink(run->GetUserOutputFileName()));
+
+  FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+  //FairLogger::GetLogger()->SetLogVerbosityLevel("VERYHIGH");
+  FairLogger::GetLogger()->SetLogVerbosityLevel("MEDIUM");
+
+  // -----   Local selection variables  -------------------------------------------
+
+  Int_t iRef    = iSel % 1000;
+  Int_t iDut    = (iSel - iRef) / 1000;
+  Int_t iDutRpc = iDut % 10;
+  iDut          = (iDut - iDutRpc) / 10;
+  Int_t iDutSm  = iDut % 10;
+  iDut          = (iDut - iDutSm) / 10;
+  Int_t iRefRpc = iRef % 10;
+  iRef          = (iRef - iRefRpc) / 10;
+  Int_t iRefSm  = iRef % 10;
+  iRef          = (iRef - iRefSm) / 10;
+
+  Int_t iSel2in  = iSel2;
+  Int_t iSel2Rpc = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Rpc) / 10;
+  Int_t iSel2Sm  = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Sm) / 10;
+
+
+  Int_t calMode = 93;
+  Int_t calSel  = 1;
+  Bool_t bOut   = kFALSE;
+
+  CbmTofEventClusterizer* tofClust =
+    new CbmTofEventClusterizer("TOF Event Clusterizer", iVerbose, bOut);
+  Int_t calSelRead = calSel;
+  if (calSel < 0) calSelRead = 0;
+  TString cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                        cFileId.Data(),
+                        iCalSet,
+                        calMode,
+                        calSelRead);
+  if (cCalId != "XXX")
+    cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                  cCalId.Data(),
+                  iCalSet,
+                  calMode,
+                  calSelRead);
+  tofClust->SetCalParFileName(cFname);
+  TString cOutFname =
+    Form("tofClust_%s_set%09d.hst.root", cFileId.Data(), iCalSet);
+  tofClust->SetOutHstFileName(cOutFname);
+
+  // =========================================================================
+  // ===                       Tracking                                    ===
+  // =========================================================================
+
+  CbmTofTrackFinder* tofTrackFinder = new CbmTofTrackFinderNN();
+  tofTrackFinder->SetMaxTofTimeDifference(0.2);  // in ns/cm
+  Int_t TrackerPar = 0;
+  switch (TrackerPar) {
+    case 0:                           // for full mTof setup
+      tofTrackFinder->SetTxLIM(0.3);  // max slope dx/dz
+      tofTrackFinder->SetTyLIM(0.3);  // max dev from mean slope dy/dz
+      tofTrackFinder->SetTxMean(0.);  // mean slope dy/dz
+      tofTrackFinder->SetTyMean(0.);  // mean slope dy/dz
+      break;
+    case 1:                             // for double stack test counters
+      tofTrackFinder->SetTxMean(0.21);  // mean slope dx/dz
+      tofTrackFinder->SetTyMean(0.18);  // mean slope dy/dz
+      tofTrackFinder->SetTxLIM(0.15);   // max slope dx/dz
+      tofTrackFinder->SetTyLIM(0.18);   // max dev from mean slope dy/dz
+      break;
+    case 2:                             // for triple stack 012 counter evaluation
+      tofTrackFinder->SetTxMean(0.0);  // mean slope dx/dz
+      tofTrackFinder->SetTyMean(-0.03);   // mean slope dy/dz
+      tofTrackFinder->SetTxLIM(0.06);  // max slope dx/dz
+      tofTrackFinder->SetTyLIM(0.03);  // max dev from mean slope dy/dz
+      break;      
+  }
+
+  CbmTofTrackFitter* tofTrackFitter = new CbmTofTrackFitterKF(0, 211);
+  TFitter* MyFit                    = new TFitter(1);  // initialize Minuit
+  tofTrackFinder->SetFitter(tofTrackFitter);
+  CbmTofFindTracks* tofFindTracks = new CbmTofFindTracks("TOF Track Finder");
+  tofFindTracks->UseFinder(tofTrackFinder);
+  tofFindTracks->UseFitter(tofTrackFitter);
+  tofFindTracks->SetCalOpt(
+    iCalOpt);  // 1 - update offsets, 2 - update walk, 0 - bypass
+  tofFindTracks->SetCorMode(iGenCor);  // valid options: 0,1,2,3,4,5,6, 10 - 19
+  //tofFindTracks->SetTtTarg(0.041);     // Mar19, Run 159
+  tofFindTracks->SetTtTarg(0.0605);  // target value for Mar2020 triple stack -> betapeak ~ 0.95
+  //tofFindTracks->SetTtTarg(0.062);   // target value for Mar2020 triple stack -> betapeak ~ 0.95
+  //tofFindTracks->SetTtTarg(0.058);   // target value for Mar2020 double stack
+  //tofFindTracks->SetTtTarg(0.051);   // target value Nov2019
+  //tofFindTracks->SetTtTarg(0.035);   // target value for inverse velocity, > 0.033 ns/cm!
+  tofFindTracks->SetCalParFileName(
+    cTrkFile);                             // Tracker parameter value file name
+  tofFindTracks->SetBeamCounter(5, 0, 0);  // default beam counter
+  tofFindTracks->SetR0Lim(0.);
+  tofFindTracks->SetStationMaxHMul(
+    30);  // Max Hit Multiplicity in any used station
+
+  tofFindTracks->SetT0MAX(dScalFac);  // in ns
+  tofFindTracks->SetSIGT(0.08);       // default in ns
+  tofFindTracks->SetSIGX(0.45);       // local-y default in cm
+  tofFindTracks->SetSIGY(0.3);        // local-x default in cm
+  tofFindTracks->SetSIGZ(10.5);       // default in cm
+  tofFindTracks->SetUseSigCalib(
+    bUseSigCalib);  // ignore resolutions in CalPar file
+  tofTrackFinder->SetSIGLIM(dChi2Lim2
+                            * 2.);     // matching window in multiples of chi2
+  tofTrackFinder->SetSIGLIMMOD(1.5);   // search window modifier for last hit (DUT)
+  tofTrackFinder->SetChiMaxAccept(dChi2Lim2);  // max tracklet chi2
+
+  Int_t iMinNofHits   = -1;
+  Int_t iNStations    = 0;
+  Int_t iNReqStations = 3;
+
+  switch (iTrackingSetup) {
+    case 0:  // bypass mode
+      iMinNofHits = -1;
+      iNStations  = 1;
+      tofFindTracks->SetStation(0, 5, 0, 0);  // Diamond
+      break;
+
+    case 1:  // for calibration mode of full setup
+    {
+      Double_t dTsig = dScalFac * 0.03;
+      tofFindTracks->SetSIGT(dTsig);  // allow for variable deviations in ns
+    }
+      iMinNofHits   = 3;
+      iNStations    = 39;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 2, 2);
+      tofFindTracks->SetStation(2, 0, 1, 2);
+      tofFindTracks->SetStation(3, 0, 0, 2);
+      tofFindTracks->SetStation(4, 0, 2, 1);
+      tofFindTracks->SetStation(5, 0, 1, 1);
+      tofFindTracks->SetStation(6, 0, 0, 1);
+      tofFindTracks->SetStation(7, 0, 2, 3);
+      tofFindTracks->SetStation(8, 0, 1, 3);
+      tofFindTracks->SetStation(9, 0, 0, 3);
+      tofFindTracks->SetStation(10, 0, 2, 0);
+      tofFindTracks->SetStation(11, 0, 1, 0);
+      tofFindTracks->SetStation(12, 0, 0, 0);
+      tofFindTracks->SetStation(13, 0, 2, 4);
+      tofFindTracks->SetStation(14, 0, 1, 4);
+      tofFindTracks->SetStation(15, 0, 0, 4);
+      tofFindTracks->SetStation(16, 0, 4, 0);
+      tofFindTracks->SetStation(17, 0, 3, 0);
+      tofFindTracks->SetStation(18, 0, 4, 1);
+      tofFindTracks->SetStation(19, 0, 3, 1);
+      tofFindTracks->SetStation(20, 0, 4, 2);
+      tofFindTracks->SetStation(21, 0, 3, 2);
+      tofFindTracks->SetStation(22, 0, 4, 3);
+      tofFindTracks->SetStation(23, 0, 3, 3);
+      tofFindTracks->SetStation(24, 0, 4, 4);
+      tofFindTracks->SetStation(25, 0, 3, 4);
+      tofFindTracks->SetStation(26, 9, 0, 0);
+      tofFindTracks->SetStation(27, 9, 0, 1);
+      tofFindTracks->SetStation(28, 7, 0, 0);
+      tofFindTracks->SetStation(29, 6, 0, 0);
+      tofFindTracks->SetStation(30, 6, 0, 1);
+      tofFindTracks->SetStation(31, 8, 0, 0);
+      tofFindTracks->SetStation(32, 8, 0, 1);
+      tofFindTracks->SetStation(33, 8, 0, 2);
+      tofFindTracks->SetStation(34, 8, 0, 3);
+      tofFindTracks->SetStation(35, 8, 0, 4);
+      tofFindTracks->SetStation(36, 8, 0, 5);
+      tofFindTracks->SetStation(37, 8, 0, 6);
+      tofFindTracks->SetStation(38, 8, 0, 7);
+      break;
+
+    case 2:
+      iMinNofHits   = 3;
+      iNStations    = 14;
+      iNReqStations = 5;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 4, 1);
+      tofFindTracks->SetStation(2, 0, 3, 1);
+      tofFindTracks->SetStation(3, 0, 4, 0);
+      tofFindTracks->SetStation(4, 0, 3, 0);
+      tofFindTracks->SetStation(5, 0, 4, 2);
+      tofFindTracks->SetStation(6, 0, 3, 2);
+      tofFindTracks->SetStation(7, 0, 4, 3);
+      tofFindTracks->SetStation(8, 0, 3, 3);
+      tofFindTracks->SetStation(9, 0, 4, 4);
+      tofFindTracks->SetStation(10, 0, 3, 4);
+      tofFindTracks->SetStation(11, 9, 0, 0);
+      tofFindTracks->SetStation(12, 9, 0, 1);
+      tofFindTracks->SetStation(13, 7, 0, 0);
+      break;
+
+    case 3:
+      iMinNofHits   = 3;
+      iNStations    = 16;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 2, 2);
+      tofFindTracks->SetStation(2, 0, 1, 2);
+      tofFindTracks->SetStation(3, 0, 0, 2);
+
+      tofFindTracks->SetStation(4, 0, 2, 1);
+      tofFindTracks->SetStation(5, 0, 1, 1);
+      tofFindTracks->SetStation(6, 0, 0, 1);
+
+      tofFindTracks->SetStation(7, 0, 2, 3);
+      tofFindTracks->SetStation(8, 0, 1, 3);
+      tofFindTracks->SetStation(9, 0, 0, 3);
+
+      tofFindTracks->SetStation(10, 0, 2, 0);
+      tofFindTracks->SetStation(11, 0, 1, 0);
+      tofFindTracks->SetStation(12, 0, 0, 0);
+
+      tofFindTracks->SetStation(13, 0, 2, 4);
+      tofFindTracks->SetStation(14, 0, 1, 4);
+      tofFindTracks->SetStation(15, 0, 0, 4);
+
+      /*
+     tofFindTracks->SetStation(16, 0, 3, 2);         
+     tofFindTracks->SetStation(17, 0, 4, 2);  
+     tofFindTracks->SetStation(18, 0, 3, 1);         
+     tofFindTracks->SetStation(19, 0, 4, 1);
+     tofFindTracks->SetStation(20, 0, 3, 3);         
+     tofFindTracks->SetStation(21, 0, 4, 3);
+     tofFindTracks->SetStation(22, 0, 3, 0);         
+     tofFindTracks->SetStation(23, 0, 4, 0);
+     tofFindTracks->SetStation(24, 0, 3, 4);         
+     tofFindTracks->SetStation(25, 0, 4, 4); 
+     */
+      break;
+
+    case 4:
+      iMinNofHits   = 3;
+      iNStations    = 4;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 2, 2);
+      tofFindTracks->SetStation(2, 0, 0, 2);
+      tofFindTracks->SetStation(3, 0, 1, 2);
+      break;
+      
+    case 14:    // for eval 012
+      iMinNofHits   = 3;
+      iNStations    = 4;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(3, 5, 0, 0);
+      tofFindTracks->SetStation(0, 0, 2, 2);
+      tofFindTracks->SetStation(1, 0, 0, 2);
+      tofFindTracks->SetStation(2, 0, 1, 2);
+      break;
+      
+    case 24:    // for geometry tuning 012
+      iMinNofHits   = 3;
+      iNStations    = 3;
+      iNReqStations = 3;
+      tofFindTracks->SetStation(0, 0, 2, 2);
+      tofFindTracks->SetStation(1, 0, 0, 2);
+      tofFindTracks->SetStation(2, 0, 1, 2);
+      break;
+
+
+    case 5:  // for calibration of 2-stack and add-on counters (STAR2, CERN)
+      iMinNofHits   = 5;
+      iNStations    = 6;
+      iNReqStations = 6;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 4, 1);
+      tofFindTracks->SetStation(2, 0, 3, 1);
+      tofFindTracks->SetStation(3, 9, 0, 0);
+      tofFindTracks->SetStation(4, 9, 0, 1);
+      tofFindTracks->SetStation(5, 7, 0, 0);
+      break;
+
+    case 6:
+      iMinNofHits   = 3;
+      iNStations    = 4;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 4, 1);
+      tofFindTracks->SetStation(2, 0, 3, 1);
+      tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+      //    tofFindTracks->SetStation(3, 9, 0, 0);
+      //    tofFindTracks->SetStation(3, 9, 0, 1);
+      //    tofFindTracks->SetStation(3, 7, 0, 0);
+      break;
+
+    case 7:  // for calibration of 2-stack and add-on counters (BUC)
+      iMinNofHits   = 3;
+      iNStations    = 5;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 4, 3);
+      tofFindTracks->SetStation(2, 0, 3, 3);
+      tofFindTracks->SetStation(3, 6, 0, 0);
+      tofFindTracks->SetStation(4, 6, 0, 1);
+      break;
+
+    case 8:  // evaluation of add-on counters (BUC)
+      iMinNofHits   = 3;
+      iNStations    = 4;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 4, 3);
+      tofFindTracks->SetStation(2, 0, 3, 3);
+      tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+      break;
+
+    case 10:
+      iMinNofHits   = 3;
+      iNStations    = 4;
+      iNReqStations = 4;
+      tofFindTracks->SetStation(0, 5, 0, 0);
+      tofFindTracks->SetStation(1, 0, 1, 2);
+      tofFindTracks->SetStation(2, 0, 0, 2);
+      tofFindTracks->SetStation(3, 0, 2, 2);
+      break;
+
+    default:
+      cout << "Tracking setup " << iTrackingSetup << " not implemented "
+           << endl;
+      return;
+      ;
+  }
+  tofFindTracks->SetMinNofHits(iMinNofHits);
+  tofFindTracks->SetNStations(iNStations);
+  tofFindTracks->SetNReqStations(iNReqStations);
+  tofFindTracks->PrintSetup();
+  run->AddTask(tofFindTracks);
+
+  // =========================================================================
+  // ===                       Analysis                                    ===
+  // =========================================================================
+  CbmTofAnaTestbeam* tofAnaTestbeam =
+    new CbmTofAnaTestbeam("TOF TestBeam Analysis", iVerbose);
+  tofAnaTestbeam->SetCorMode(iAnaCor);  // 1 - DTD4, 2 - X4, 3 - Y4, 4 - Texp
+  tofAnaTestbeam->SetHitDistMin(30.);   // initialization
+  tofAnaTestbeam->SetEnableMatchPosScaling(kTRUE);
+  tofAnaTestbeam->SetSpillDuration(3.);
+  if (iMc == 1) {
+    tofAnaTestbeam->SetSpillDuration(0.);
+    tofAnaTestbeam->SetSpillBreak(0.);
+  }
+  //CbmTofAnaTestbeam defaults
+  tofAnaTestbeam->SetR0LimFit(
+    0.);  // limit distance of fitted track to nominal vertex
+  tofAnaTestbeam->SetDXMean(0.);
+  tofAnaTestbeam->SetDYMean(0.);
+  tofAnaTestbeam->SetDTMean(0.);  // in ns
+  tofAnaTestbeam->SetDXWidth(0.5);
+  tofAnaTestbeam->SetDYWidth(1.0);
+  tofAnaTestbeam->SetDTWidth(0.1);  // in ns
+  tofAnaTestbeam->SetCalParFileName(cAnaFile);
+  Double_t dScalFacA = 0.9;  // dScalFac is used for tracking
+  tofAnaTestbeam->SetPosY4Sel(
+    0.5 * dScalFacA);  // Y Position selection in fraction of strip length
+  tofAnaTestbeam->SetDTDia(0.);    // Time difference to additional diamond
+  tofAnaTestbeam->SetMul0Max(20);  // Max Multiplicity in dut
+  tofAnaTestbeam->SetMul4Max(30);  // Max Multiplicity in Ref - RPC
+  tofAnaTestbeam->SetMulDMax(3);   // Max Multiplicity in Diamond / BeamRef
+  tofAnaTestbeam->SetTOffD4(14.);  // initialization
+  tofAnaTestbeam->SetDTD4MAX(
+    6.);  // initialization of Max time difference Ref - BRef
+
+  //tofAnaTestbeam->SetTShift(-28000.);// initialization
+  tofAnaTestbeam->SetPosYS2Sel(
+    0.55);  // Y Position selection in fraction of strip length
+  tofAnaTestbeam->SetChS2Sel(0.);     // Center of channel selection window
+  tofAnaTestbeam->SetDChS2Sel(100.);  // Width  of channel selection window
+  tofAnaTestbeam->SetSel2TOff(0.);    // Shift Sel2 time peak to 0
+  tofAnaTestbeam->SetChi2Lim(5.);     // initialization of Chi2 selection limit
+  tofAnaTestbeam->SetChi2Lim2(
+    3.);  // initialization of Chi2 selection limit for Mref-Sel2 pair
+  tofAnaTestbeam->SetDutDX(
+    15.);  // limit inspection of tracklets to selected region
+  tofAnaTestbeam->SetDutDY(
+    15.);  // limit inspection of tracklets to selected region
+  tofAnaTestbeam->SetSIGLIM(3.);  // max matching chi2
+  tofAnaTestbeam->SetSIGT(0.08);  // in ns
+  tofAnaTestbeam->SetSIGX(0.3);   // in cm
+  tofAnaTestbeam->SetSIGY(0.6);   // in cm
+
+  Int_t iRSel    = 500;
+  Int_t iRSelTyp = 5;
+  Int_t iRSelSm  = 0;
+  Int_t iRSelRpc = 0;
+  Int_t iRSelin  = iRSel;
+
+  tofAnaTestbeam->SetBeamRefSmType(iRSelTyp);  // common reaction reference
+  tofAnaTestbeam->SetBeamRefSmId(iRSelSm);
+  tofAnaTestbeam->SetBeamRefRpc(iRSelRpc);
+
+  if (iSel2 >= -1) {
+    tofAnaTestbeam->SetMrpcSel2(
+      iSel2);  // initialization of second selector Mrpc Type
+    tofAnaTestbeam->SetMrpcSel2Sm(
+      iSel2Sm);  // initialization of second selector Mrpc SmId
+    tofAnaTestbeam->SetMrpcSel2Rpc(
+      iSel2Rpc);  // initialization of second selector Mrpc RpcId
+  }
+
+  cout << "AnaTestbeam init for Dut " << iDut << iDutSm << iDutRpc << ", Ref "
+       << iRef << iRefSm << iRefRpc << endl;
+
+  tofAnaTestbeam->SetDut(iDut);            // Device under test
+  tofAnaTestbeam->SetDutSm(iDutSm);        // Device under test
+  tofAnaTestbeam->SetDutRpc(iDutRpc);      // Device under test
+  tofAnaTestbeam->SetMrpcRef(iRef);        // Reference RPC
+  tofAnaTestbeam->SetMrpcRefSm(iRefSm);    // Reference RPC
+  tofAnaTestbeam->SetMrpcRefRpc(iRefRpc);  // Reference RPC
+
+  cout << "dispatch iSel = " << iSel << ", iSel2in = " << iSel2in
+       << ", iRSelin = " << iRSelin << ", iRSel = " << iRSel << endl;
+  if (1) {
+    switch (iSel) {
+
+      case 10:
+        switch (iRSelin) {
+          case 500:
+            tofAnaTestbeam->SetTShift(2.5);  // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(18.);  // Shift DTD4 to physical value
+            switch (iSel2in) {
+              case 20:
+                tofAnaTestbeam->SetSel2TOff(0.);  // Shift Sel2 time peak to 0
+                break;
+              default:;
+            }
+            break;
+          default:;
+        }
+        break;
+
+      case 700040:
+      case 900040:
+      case 901040:
+        switch (iRSelin) {
+          case 500:
+            tofAnaTestbeam->SetTShift(0.3);  // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(18.);  // Shift DTD4 to physical value
+
+            switch (iSel2in) {
+              case 30:
+                tofAnaTestbeam->SetSel2TOff(-0.3);  // Shift Sel2 time peak to 0
+                break;
+              case 31:
+                tofAnaTestbeam->SetSel2TOff(
+                  -0.41);  // Shift Sel2 time peak to 0
+                break;
+
+              default:;
+            }
+            break;
+          default:;
+        }
+        break;
+
+      case 700041:
+      case 900041:
+      case 901041:
+        switch (iRSelin) {
+          case 500:
+            tofAnaTestbeam->SetTShift(0.8);  // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(11.);  // Shift DTD4 to physical value
+
+            switch (iSel2in) {
+              case 30:
+                tofAnaTestbeam->SetSel2TOff(-0.3);  // Shift Sel2 time peak to 0
+                break;
+              case 31:
+                tofAnaTestbeam->SetSel2TOff(0.);  // Shift Sel2 time peak to 0
+                break;
+              case 600:
+                tofAnaTestbeam->SetSel2TOff(-0.2);  // Shift Sel2 time peak to 0
+                break;
+
+              default:;
+            }
+            break;
+          default:;
+        }
+        break;
+
+      case 600043:
+      case 601043:
+        switch (iRSelin) {
+          case 500:
+            tofAnaTestbeam->SetTShift(5.3);  // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(11.);  // Shift DTD4 to physical value
+
+            switch (iSel2in) {
+              case 33:
+                tofAnaTestbeam->SetSel2TOff(
+                  -0.55);  // Shift Sel2 time peak to 0
+                break;
+
+              default:;
+            }
+            break;
+          default:;
+        }
+        break;
+
+      default:
+        cout << "Better to define analysis setup! Running with default offset "
+                "parameter... "
+             << endl;
+        // return;
+    }  // end of different subsets
+
+    cout << " Initialize TSHIFT to " << tofAnaTestbeam->GetTShift() << endl;
+    run->AddTask(tofAnaTestbeam);
+  }
+
+  // -----  Parameter database   --------------------------------------------
+  FairRuntimeDb* rtdb       = run->GetRuntimeDb();
+  Bool_t kParameterMerged   = kTRUE;
+  FairParRootFileIo* parIo2 = new FairParRootFileIo(kParameterMerged);
+  parIo2->open(ParFile.Data(), "UPDATE");
+  parIo2->print();
+  rtdb->setFirstInput(parIo2);
+
+  FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
+  parIo1->open(parFileList, "in");
+  parIo1->print();
+  rtdb->setSecondInput(parIo1);
+  rtdb->print();
+  rtdb->printParamContexts();
+
+  //  FairParRootFileIo* parInput1 = new FairParRootFileIo();
+  //  parInput1->open(ParFile.Data());
+  //  rtdb->setFirstInput(parInput1);
+
+  // -----   Intialise and run   --------------------------------------------
+  run->Init();
+  cout << "Starting run" << endl;
+  run->Run(0, nEvents);
+  //run->Run(nEvents-1, nEvents); //debugging single events for memory leak
+  // ------------------------------------------------------------------------
+
+  // default displays, plot results
+  /*
+  TString Display_Status = "pl_over_Mat04D4best.C";
+  TString Display_Funct;
+  if (iGenCor<0) {
+    Display_Funct = "pl_over_Mat04D4best(1)";
+  }else{
+    Display_Funct = "pl_over_Mat04D4best(0)";
+  }
+  gROOT->LoadMacro(Display_Status);
+
+  cout << "Exec "<< Display_Funct.Data()<< endl;
+  gInterpreter->ProcessLine(Display_Funct);
+  */
+  gROOT->LoadMacro("save_hst.C");
+  gROOT->LoadMacro("pl_over_MatD4sel.C");
+  gROOT->LoadMacro("pl_eff_XY.C");
+  gROOT->LoadMacro("pl_over_trk.C");
+  gROOT->LoadMacro("pl_calib_trk.C");
+  gROOT->LoadMacro("pl_XY_trk.C");
+  gROOT->LoadMacro("pl_vert_trk.C");
+  gROOT->LoadMacro("pl_pull_trk.C");
+  gROOT->LoadMacro("pl_all_Track2D.C");
+  gROOT->LoadMacro("pl_TIS.C");
+  gROOT->LoadMacro("pl_TIR.C");
+  gROOT->LoadMacro("pl_Eff_XY.C");
+  gROOT->LoadMacro("pl_Eff_DTLH.C");
+  gROOT->LoadMacro("pl_Eff_TIS.C");
+  gROOT->LoadMacro("pl_Dut_Res.C");
+
+  TString SaveToHstFile = "save_hst(\"" + cHstFile + "\")";
+  gInterpreter->ProcessLine(SaveToHstFile);
+  
+  return;
+
+  //gInterpreter->ProcessLine("pl_over_MatD4sel()");
+  //gInterpreter->ProcessLine("pl_TIS()");
+  //gInterpreter->ProcessLine("pl_TIR()");
+  //gInterpreter->ProcessLine("pl_eff_XY()");
+  gInterpreter->ProcessLine("pl_calib_trk()");
+  gInterpreter->ProcessLine("pl_vert_trk()");
+
+  gInterpreter->ProcessLine("pl_all_Track2D(1)");
+  gInterpreter->ProcessLine("pl_all_Track2D(2)");
+  gInterpreter->ProcessLine("pl_all_Track2D(4)");
+
+  TString over_trk = "pl_over_trk(" + (TString)(Form("%d", iNStations)) + ")";
+  gInterpreter->ProcessLine(over_trk);
+
+  TString XY_trk = "pl_XY_trk(" + (TString)(Form("%d", iNStations)) + ")";
+  gInterpreter->ProcessLine(XY_trk);
+
+  TString Pull0 = (TString)(Form("pl_pull_trk(%d,%d,1)", iNStations, 0));
+  gInterpreter->ProcessLine(Pull0);
+  TString Pull1 = (TString)(Form("pl_pull_trk(%d,%d,1)", iNStations, 1));
+  gInterpreter->ProcessLine(Pull1);
+  TString Pull3 = (TString)(Form("pl_pull_trk(%d,%d,1)", iNStations, 3));
+  gInterpreter->ProcessLine(Pull3);
+  TString Pull4 = (TString)(Form("pl_pull_trk(%d,%d,1)", iNStations, 4));
+  gInterpreter->ProcessLine(Pull4);
+}
diff --git a/macro/beamtime/mcbm2018/dis_digi.C b/macro/beamtime/mcbm2018/dis_digi.C
new file mode 100644
index 0000000000000000000000000000000000000000..272f454ffa88c9f95a0d172cac2ccd8647a9f9f5
--- /dev/null
+++ b/macro/beamtime/mcbm2018/dis_digi.C
@@ -0,0 +1,860 @@
+void dis_digi(Int_t nEvents        = 100,
+              Int_t calMode        = 93,
+              Int_t calSel         = 1,
+              Int_t calSm          = 0,
+              Int_t RefSel         = 1,
+              TString cFileId      = "159.50.4.1",
+              Int_t iCalSet        = 10500,
+              Bool_t bOut          = 0,
+              Int_t iSel2          = 20,
+              Double_t dDeadtime   = 50,
+              Int_t iGenCor        = 1,
+              Int_t iTrackingSetup = 1,
+              Double_t dScalFac    = 5.,
+              Double_t dChi2Lim2   = 10.,
+              TString cCalId       = "XXX",
+              Bool_t bUseSigCalib  = kFALSE,
+              Int_t iCalOpt        = 1) {
+  
+  Int_t iVerbose = 1;
+  if (cCalId == "") cCalId = cFileId;
+  TString FId = cFileId;
+  TString cRun(FId(0, 3));
+  Int_t iRun = cRun.Atoi();
+  cout << "dis_digi for Run " << iRun << endl;
+  
+  //Specify log level (INFO, DEBUG, DEBUG1, ...)
+  //TString logLevel = "FATAL";
+  //TString logLevel = "ERROR";
+  //TString logLevel = "INFO";
+  TString logLevel = "DEBUG";
+  //TString logLevel = "DEBUG1";
+  //TString logLevel = "DEBUG2";
+  //TString logLevel = "DEBUG3";
+  FairLogger::GetLogger();
+  gLogger->SetLogScreenLevel(logLevel);
+  //gLogger->SetLogScreenLevel("DEBUG");
+  gLogger->SetLogVerbosityLevel("MEDIUM");
+
+  TString workDir = gSystem->Getenv("VMCWORKDIR");
+  /*
+   TString workDir    = (TString)gInterpreter->ProcessLine(".! pwd");
+   cout << "workdir = "<< workDir.Data() << endl;
+   return;
+  */
+  //   TString paramDir   = workDir + "/macro/beamtime/mcbm2018/";
+  TString paramDir   = "./";
+  TString ParFile    = paramDir + "data/" + cFileId + ".params.root";
+  TString InputFile  = paramDir + "data/" + cFileId + ".root";
+  TString OutputFile = paramDir + "data/disdigi_" + cFileId
+                       + Form("_%09d%03d", iCalSet, iSel2) + ".out.root";
+
+  TString cTrkFile = Form("%s_tofFindTracks.hst.root", cFileId.Data());
+
+  TList* parFileList = new TList();
+
+
+  TString shcmd = "rm -v " + ParFile;
+  gSystem->Exec(shcmd.Data());
+  
+  Int_t iGeo = 0;  //iMc;
+  if (iGeo == 0) {
+    TString TofGeo = "v18m_mcbm";
+
+    cout << "Geometry version " << TofGeo << endl;
+    /*
+    TObjString* tofDigiFile = new TObjString(
+      workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par");  // TOF digi file
+    parFileList->Add(tofDigiFile);
+    */
+    TObjString* tofDigiBdfFile =
+      new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digibdf.par");
+    parFileList->Add(tofDigiBdfFile);
+
+    TString geoDir  = gSystem->Getenv("VMCWORKDIR");
+    TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+    TFile* fgeo     = new TFile(geoFile);
+    TGeoManager* geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+    if (NULL == geoMan) {
+      cout << "<E> FAIRGeom not found in geoFile" << endl;
+      return;
+    }
+
+    if (0) {
+      TGeoVolume* master = geoMan->GetTopVolume();
+      master->SetVisContainers(1);
+      master->Draw("ogl");
+    }
+  }
+  
+  // Local steering variables
+  Int_t iBRef    = iCalSet % 1000;
+  Int_t iSet     = (iCalSet - iBRef) / 1000;
+  Int_t iRSel    = 0;
+  Int_t iRSelTyp = 0;
+  Int_t iRSelSm  = 0;
+  Int_t iRSelRpc = 0;
+  if (iSel2 == 0) {
+    iRSel = iBRef;  // use diamond
+    iSel2 = iBRef;
+  } else {
+    if (iSel2 < 0) iSel2 = -iSel2;
+    iRSel = iSel2;
+  }
+
+  iRSelRpc = iRSel % 10;
+  iRSelTyp = (iRSel - iRSelRpc) / 10;
+  iRSelSm  = iRSelTyp % 10;
+  iRSelTyp = (iRSelTyp - iRSelSm) / 10;
+
+  Int_t iSel2in  = iSel2;
+  Int_t iSel2Rpc = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Rpc) / 10;
+  Int_t iSel2Sm  = iSel2 % 10;
+  iSel2          = (iSel2 - iSel2Sm) / 10;
+
+  Int_t iRef    = iSet % 1000;
+  Int_t iDut    = (iSet - iRef) / 1000;
+  Int_t iDutRpc = iDut % 10;
+  iDut          = (iDut - iDutRpc) / 10;
+  Int_t iDutSm  = iDut % 10;
+  iDut          = (iDut - iDutSm) / 10;
+
+  Int_t iRefRpc = iRef % 10;
+  iRef          = (iRef - iRefRpc) / 10;
+  Int_t iRefSm  = iRef % 10;
+  iRef          = (iRef - iRefSm) / 10;
+
+  // -----   Reconstruction run   -------------------------------------------
+  FairRunAna* run = new FairRunAna();
+  run->SetInputFile(InputFile.Data());
+  //run->SetOutputFile(OutputFile);
+  run->SetUserOutputFileName(OutputFile.Data());
+  run->SetSink(new FairRootFileSink(run->GetUserOutputFileName()));
+
+  CbmTofEventClusterizer* tofClust =
+    new CbmTofEventClusterizer("TOF Event Clusterizer", iVerbose, bOut);
+
+  tofClust->SetCalMode(calMode);
+  tofClust->SetCalSel(calSel);
+  tofClust->SetCaldXdYMax(50.);  // geometrical matching window in cm
+  tofClust->SetCalCluMulMax(
+    4.);  // Max Counter Cluster Multiplicity for filling calib histos
+  tofClust->SetCalRpc(calSm);   // select detector for calibration update
+  tofClust->SetTRefId(RefSel);  // reference trigger for offset calculation
+  tofClust->SetTotMax(20.);     // Tot upper limit for walk corection
+  tofClust->SetTotMin(
+    0.01);  //(12000.);  // Tot lower limit for walk correction
+  tofClust->SetTotPreRange(
+    2.);  // effective lower Tot limit  in ns from peak position
+  tofClust->SetTotMean(2.);       // Tot calibration target value in ns
+  tofClust->SetMaxTimeDist(1.0);  // default cluster range in ns
+  //tofClust->SetMaxTimeDist(0.);       //Deb// default cluster range in ns
+  tofClust->SetDelTofMax(
+    60.);  // acceptance range for cluster correlation in cm (!)
+  tofClust->SetSel2MulMax(4);  // limit Multiplicity in 2nd selector
+  tofClust->SetChannelDeadtime(dDeadtime);  // artificial deadtime in ns
+  tofClust->SetEnableAvWalk(kTRUE);
+  tofClust->SetYFitMin(1.E4);
+  //tofClust->SetTimePeriod(6.25);       // ignore coarse time
+  //tofClust->SetCorMode(2);              // correct missing hits
+
+  Int_t calSelRead = calSel;
+  if (calSel < 0) calSelRead = 0;
+  TString cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                        cFileId.Data(),
+                        iCalSet,
+                        calMode,
+                        calSelRead);
+  tofClust->SetCalParFileName(cFname);
+  TString cOutFname =
+    Form("tofClust_%s_set%09d.hst.root", cFileId.Data(), iCalSet);
+  tofClust->SetOutHstFileName(cOutFname);
+
+  TString cAnaFile =
+    Form("%s_%09d%03d_tofAna.hst.root", cFileId.Data(), iCalSet, iSel2);
+
+  switch (calMode) {
+    case -1:                      // initial calibration
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      //tofClust->SetTotMin(1.);
+      tofClust->SetTRefDifMax(50000.);  // in ns
+      tofClust->PosYMaxScal(10000.);    // in % of length
+      tofClust->SetMaxTimeDist(0.);     // no cluster building
+      //tofClust->SetTimePeriod(0.);       // inspect coarse time
+      break;
+    case 0:                       // initial calibration
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      //tofClust->SetTotMin(1.);
+      tofClust->SetTRefDifMax(50.);  // in ns
+      tofClust->PosYMaxScal(10.);    // in % of length
+      tofClust->SetMaxTimeDist(0.);  // no cluster building
+      //tofClust->SetTimePeriod(0.);       // inspect coarse time
+      break;
+    case 1:                       // save offsets, update walks, for diamonds
+      tofClust->SetTotMax(256.);  // range in bin number
+      tofClust->SetTotPreRange(256.);
+      tofClust->SetTRefDifMax(6.25);  // in ns
+      //tofClust->SetTimePeriod(6.25);       // inspect coarse time
+      tofClust->PosYMaxScal(10.);  // in % of length
+      break;
+    case 11:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(3.0);   // in % of length
+      break;
+    case 21:
+      tofClust->SetTRefDifMax(2.5);  // in ns
+      tofClust->PosYMaxScal(2.0);    // in % of length
+      break;
+    case 31:
+      tofClust->SetTRefDifMax(2.);  // in ns
+      tofClust->PosYMaxScal(1.5);   // in % of length
+      break;
+    case 41:
+      tofClust->SetTRefDifMax(1.);  // in ns
+      tofClust->PosYMaxScal(0.8);   // in % of length
+      break;
+    case 51:
+      tofClust->SetTRefDifMax(0.7);  // in ns
+      tofClust->PosYMaxScal(0.7);    // in % of length
+      break;
+    case 61:
+      tofClust->SetTRefDifMax(0.5);  // in ns
+      tofClust->PosYMaxScal(0.7);    // in % of length
+      break;
+    case 71:
+      tofClust->SetTRefDifMax(0.4);  // in ns
+      tofClust->PosYMaxScal(0.6);    // in % of length
+      break;
+
+    case 2:                           // time difference calibration
+      tofClust->SetTRefDifMax(300.);  // in ns
+      tofClust->PosYMaxScal(1000.);   //in % of length
+      break;
+
+    case 3:                           // time offsets
+      tofClust->SetTRefDifMax(200.);  // in ns
+      tofClust->PosYMaxScal(100.);    //in % of length
+      tofClust->SetMaxTimeDist(0.);   // no cluster building
+      break;
+    case 12:
+    case 13:
+      tofClust->SetTRefDifMax(100.);  // in ns
+      tofClust->PosYMaxScal(50.);     //in % of length
+      break;
+    case 22:
+    case 23:
+      tofClust->SetTRefDifMax(50.);  // in ns
+      tofClust->PosYMaxScal(20.);    //in % of length
+      break;
+    case 32:
+    case 33:
+      tofClust->SetTRefDifMax(25.);  // in ns
+      tofClust->PosYMaxScal(10.);    //in % of length
+      break;
+    case 42:
+    case 43:
+      tofClust->SetTRefDifMax(12.);  // in ns
+      tofClust->PosYMaxScal(5.);     //in % of length
+      break;
+    case 52:
+    case 53:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(3.);    //in % of length
+      break;
+    case 62:
+    case 63:
+      tofClust->SetTRefDifMax(3.);  // in ns
+      tofClust->PosYMaxScal(2.);    //in % of length
+      break;
+    case 72:
+    case 73:
+      tofClust->SetTRefDifMax(2.);  // in ns
+      tofClust->PosYMaxScal(1.5);   //in % of length
+      break;
+    case 82:
+    case 83:
+      tofClust->SetTRefDifMax(1.);  // in ns
+      tofClust->PosYMaxScal(1.0);   //in % of length
+      break;
+    case 92:
+    case 93:
+      tofClust->SetTRefDifMax(0.6);  // in ns
+      tofClust->PosYMaxScal(1.0);    //in % of length
+      break;
+
+    case 4:                         // velocity dependence (DelTOF)
+      tofClust->SetTRefDifMax(6.);  // in ns
+      tofClust->PosYMaxScal(1.5);   //in % of length
+      break;
+    case 14:
+      tofClust->SetTRefDifMax(5.);  // in ns
+      tofClust->PosYMaxScal(1.);    //in % of length
+      break;
+    case 24:
+      tofClust->SetTRefDifMax(3.);  // in ns
+      tofClust->PosYMaxScal(1.0);   //in % of length
+      break;
+    case 34:
+      tofClust->SetTRefDifMax(2.);  // in ns
+      tofClust->PosYMaxScal(1.0);   //in % of length
+      break;
+    case 44:
+      tofClust->SetTRefDifMax(1.);  // in ns
+      tofClust->PosYMaxScal(1.0);   //in % of length
+      break;
+    case 54:
+      tofClust->SetTRefDifMax(0.7);  // in ns
+      tofClust->PosYMaxScal(0.7);    //in % of length
+      break;
+    case 64:
+      tofClust->SetTRefDifMax(0.5);  // in ns
+      tofClust->PosYMaxScal(0.7);    //in % of length
+      break;
+    default:
+      cout << "<E> Calib mode not implemented! stop execution of script"
+           << endl;
+      return;
+  }
+
+  run->AddTask(tofClust);
+
+  // =========================================================================
+  // ===                       Tracking                                    ===
+  // =========================================================================
+
+  CbmTofTrackFinder* tofTrackFinder = new CbmTofTrackFinderNN();
+    tofTrackFinder->SetMaxTofTimeDifference(0.4);  // in ns/cm
+    Int_t TrackerPar = 0;
+    switch (TrackerPar) {
+      case 0:                           // for full mTof setup
+        tofTrackFinder->SetTxLIM(0.3);  // max slope dx/dz
+        tofTrackFinder->SetTyLIM(0.3);  // max dev from mean slope dy/dz
+        tofTrackFinder->SetTxMean(0.);  // mean slope dy/dz
+        tofTrackFinder->SetTyMean(0.);  // mean slope dy/dz
+        break;
+      case 1:                             // for double stack test counters
+        tofTrackFinder->SetTxMean(0.21);  // mean slope dy/dz
+        tofTrackFinder->SetTyMean(0.18);  // mean slope dy/dz
+        tofTrackFinder->SetTxLIM(0.15);   // max slope dx/dz
+        tofTrackFinder->SetTyLIM(0.18);   // max dev from mean slope dy/dz
+        break;
+    }
+
+    CbmTofTrackFitter* tofTrackFitter = new CbmTofTrackFitterKF(0, 211);
+    TFitter* MyFit                    = new TFitter(1);  // initialize Minuit
+    tofTrackFinder->SetFitter(tofTrackFitter);
+    CbmTofFindTracks* tofFindTracks = new CbmTofFindTracks("TOF Track Finder");
+    tofFindTracks->UseFinder(tofTrackFinder);
+    tofFindTracks->UseFitter(tofTrackFitter);
+    tofFindTracks->SetCalOpt(
+      iCalOpt);  // 1 - update offsets, 2 - update walk, 0 - bypass
+    tofFindTracks->SetCorMode(
+      iGenCor);  // valid options: 0,1,2,3,4,5,6, 10 - 19
+    tofFindTracks->SetTtTarg(
+      0.057);  // target value for inverse velocity, > 0.033 ns/cm!
+    //tofFindTracks->SetTtTarg(0.035);            // target value for inverse velocity, > 0.033 ns/cm!
+    tofFindTracks->SetCalParFileName(
+      cTrkFile);  // Tracker parameter value file name
+    tofFindTracks->SetBeamCounter(5, 0, 0);  // default beam counter
+    tofFindTracks->SetR0Lim(100.);
+
+    tofFindTracks->SetStationMaxHMul(
+      30);  // Max Hit Multiplicity in any used station
+
+    tofFindTracks->SetT0MAX(dScalFac);  // in ns
+    tofFindTracks->SetSIGT(0.08);       // default in ns
+    tofFindTracks->SetSIGX(0.3);        // default in cm
+    tofFindTracks->SetSIGY(0.6);        // default in cm
+    tofFindTracks->SetSIGZ(0.05);       // default in cm
+    tofFindTracks->SetUseSigCalib(
+      bUseSigCalib);  // ignore resolutions in CalPar file
+    tofTrackFinder->SetSIGLIM(dChi2Lim2
+                              * 2.);  // matching window in multiples of chi2
+    tofTrackFinder->SetChiMaxAccept(dChi2Lim2);  // max tracklet chi2
+
+
+    Int_t iMinNofHits   = -1;
+    Int_t iNStations    = 0;
+    Int_t iNReqStations = 3;
+
+    switch (iTrackingSetup) {
+      case 0:  // bypass mode
+        iMinNofHits = -1;
+        iNStations  = 1;
+        tofFindTracks->SetStation(0, 5, 0, 0);  // Diamond
+        break;
+
+      case 1:  // for calibration mode of full setup
+      {
+        Double_t dTsig = dScalFac * 0.03;
+        tofFindTracks->SetSIGT(dTsig);  // allow for variable deviations in ns
+      }
+        iMinNofHits   = 3;
+        iNStations    = 28;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 2, 2);
+        tofFindTracks->SetStation(2, 0, 1, 2);
+        tofFindTracks->SetStation(3, 0, 0, 2);
+        tofFindTracks->SetStation(4, 0, 2, 1);
+        tofFindTracks->SetStation(5, 0, 1, 1);
+        tofFindTracks->SetStation(6, 0, 0, 1);
+        tofFindTracks->SetStation(7, 0, 2, 3);
+        tofFindTracks->SetStation(8, 0, 1, 3);
+        tofFindTracks->SetStation(9, 0, 0, 3);
+        tofFindTracks->SetStation(10, 0, 2, 0);
+        tofFindTracks->SetStation(11, 0, 1, 0);
+        tofFindTracks->SetStation(12, 0, 0, 0);
+        tofFindTracks->SetStation(13, 0, 2, 4);
+        tofFindTracks->SetStation(14, 0, 1, 4);
+        tofFindTracks->SetStation(15, 0, 0, 4);
+        tofFindTracks->SetStation(16, 0, 4, 0);
+        tofFindTracks->SetStation(17, 0, 3, 0);
+        tofFindTracks->SetStation(18, 0, 4, 1);
+        tofFindTracks->SetStation(19, 0, 3, 1);
+        tofFindTracks->SetStation(20, 0, 4, 2);
+        tofFindTracks->SetStation(21, 0, 3, 2);
+        tofFindTracks->SetStation(22, 0, 4, 3);
+        tofFindTracks->SetStation(23, 0, 3, 3);
+        tofFindTracks->SetStation(24, 0, 4, 4);
+        tofFindTracks->SetStation(25, 0, 3, 4);
+        tofFindTracks->SetStation(26, 9, 0, 0);
+        tofFindTracks->SetStation(27, 9, 0, 1);
+        //tofFindTracks->SetStation(28, 6, 0, 0);
+        //tofFindTracks->SetStation(29, 6, 0, 1);
+        break;
+
+      case 2:
+        iMinNofHits   = 3;
+        iNStations    = 14;
+        iNReqStations = 5;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 1);
+        tofFindTracks->SetStation(2, 0, 3, 1);
+        tofFindTracks->SetStation(3, 0, 4, 0);
+        tofFindTracks->SetStation(4, 0, 3, 0);
+        tofFindTracks->SetStation(5, 0, 4, 2);
+        tofFindTracks->SetStation(6, 0, 3, 2);
+        tofFindTracks->SetStation(7, 0, 4, 3);
+        tofFindTracks->SetStation(8, 0, 3, 3);
+        tofFindTracks->SetStation(9, 0, 4, 4);
+        tofFindTracks->SetStation(10, 0, 3, 4);
+        tofFindTracks->SetStation(11, 9, 0, 0);
+        tofFindTracks->SetStation(12, 9, 0, 1);
+        tofFindTracks->SetStation(13, 7, 0, 0);
+        break;
+
+      case 3:
+        iMinNofHits   = 3;
+        iNStations    = 16;
+        iNReqStations = 4;
+        
+		tofFindTracks->SetStation(0, 0, 2, 2);
+		tofFindTracks->SetStation(1, 0, 1, 2);
+		tofFindTracks->SetStation(2, 0, 0, 2);
+
+		tofFindTracks->SetStation(3, 0, 2, 1);
+		tofFindTracks->SetStation(4, 0, 1, 1);
+		tofFindTracks->SetStation(5, 0, 0, 1);
+
+		tofFindTracks->SetStation(6, 0, 2, 3);
+		tofFindTracks->SetStation(7, 0, 1, 3);
+		tofFindTracks->SetStation(8, 0, 0, 3);
+
+		tofFindTracks->SetStation(9,  0, 2, 0);
+		tofFindTracks->SetStation(10, 0, 1, 0);
+		tofFindTracks->SetStation(11, 0, 0, 0);
+
+		tofFindTracks->SetStation(12, 0, 2, 4);
+		tofFindTracks->SetStation(13, 0, 1, 4);
+		tofFindTracks->SetStation(14, 0, 0, 4);
+
+		tofFindTracks->SetStation(15, 5, 0, 0);
+        /*
+     tofFindTracks->SetStation(16, 0, 3, 2);         
+     tofFindTracks->SetStation(17, 0, 4, 2);  
+     tofFindTracks->SetStation(18, 0, 3, 1);         
+     tofFindTracks->SetStation(19, 0, 4, 1);
+     tofFindTracks->SetStation(20, 0, 3, 3);         
+     tofFindTracks->SetStation(21, 0, 4, 3);
+     tofFindTracks->SetStation(22, 0, 3, 0);         
+     tofFindTracks->SetStation(23, 0, 4, 0);
+     tofFindTracks->SetStation(24, 0, 3, 4);         
+     tofFindTracks->SetStation(25, 0, 4, 4); 
+     */
+        break;
+
+      case 4:
+        iMinNofHits   = 3;
+        iNStations    = 4;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 1);
+        tofFindTracks->SetStation(2, 0, 3, 1);
+        tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+        break;
+
+      case 5:  // for calibration of 2-stack and add-on counters (STAR2, BUC)
+        iMinNofHits   = 3;
+        iNStations    = 5;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 1);
+        tofFindTracks->SetStation(2, 0, 3, 1);
+        tofFindTracks->SetStation(3, 9, 0, 0);
+        tofFindTracks->SetStation(4, 9, 0, 1);
+        break;
+
+      case 6:
+        iMinNofHits   = 3;
+        iNStations    = 4;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 1);
+        tofFindTracks->SetStation(2, 0, 3, 1);
+        tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+        //    tofFindTracks->SetStation(3, 9, 0, 0);
+        //    tofFindTracks->SetStation(3, 9, 0, 1);
+        //    tofFindTracks->SetStation(3, 7, 0, 0);
+        break;
+
+      case 7:  // for calibration of 2-stack and add-on counters (BUC)
+        iMinNofHits   = 4;
+        iNStations    = 5;
+        iNReqStations = 5;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 3);
+        tofFindTracks->SetStation(2, 0, 3, 3);
+        tofFindTracks->SetStation(3, 6, 0, 0);
+        tofFindTracks->SetStation(4, 6, 0, 1);
+        break;
+
+      case 8:  // evaluation of add-on counters (BUC)
+        iMinNofHits   = 3;
+        iNStations    = 4;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 3);
+        tofFindTracks->SetStation(2, 0, 3, 3);
+        tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+        break;
+
+      case 10:
+        iMinNofHits   = 3;
+        iNStations    = 4;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 1, 2);
+        tofFindTracks->SetStation(2, 0, 0, 2);
+        tofFindTracks->SetStation(3, 0, 2, 2);
+
+      default:
+        cout << "Tracking setup " << iTrackingSetup << " not implemented "
+             << endl;
+        return;
+        ;
+    }
+    tofFindTracks->SetMinNofHits(iMinNofHits);
+    tofFindTracks->SetNStations(iNStations);
+    tofFindTracks->SetNReqStations(iNReqStations);
+    tofFindTracks->PrintSetup();
+    run->AddTask(tofFindTracks);
+
+    // =========================================================================
+    // ===                       Analysis                                    ===
+    // =========================================================================
+
+  CbmTofAnaTestbeam* tofAnaTestbeam =
+    new CbmTofAnaTestbeam("TOF TestBeam Analysis", iVerbose);
+
+  //CbmTofAnaTestbeam defaults
+  tofAnaTestbeam->SetReqTrg(0);  // 0 - no selection
+  tofAnaTestbeam->SetDXMean(0.);
+  tofAnaTestbeam->SetDYMean(0.);
+  tofAnaTestbeam->SetDTMean(0.);  // in ps
+  tofAnaTestbeam->SetDXWidth(0.4);
+  tofAnaTestbeam->SetDYWidth(0.4);
+  tofAnaTestbeam->SetDTWidth(80.);  // in ps
+  tofAnaTestbeam->SetCalParFileName(cAnaFile);
+  tofAnaTestbeam->SetPosY4Sel(
+    0.5);  // Y Position selection in fraction of strip length
+  tofAnaTestbeam->SetDTDia(0.);        // Time difference to additional diamond
+  tofAnaTestbeam->SetCorMode(RefSel);  // 1 - DTD4, 2 - X4
+  tofAnaTestbeam->SetMul0Max(30);      // Max Multiplicity in dut
+  tofAnaTestbeam->SetMul4Max(30);      // Max Multiplicity in Ref - RPC
+  tofAnaTestbeam->SetMulDMax(10);      // Max Multiplicity in Diamond
+  tofAnaTestbeam->SetHitDistMin(30.);  // initialization
+  tofAnaTestbeam->SetEnableMatchPosScaling(kTRUE);
+
+  tofAnaTestbeam->SetPosYS2Sel(
+    0.5);  // Y Position selection in fraction of strip length
+  tofAnaTestbeam->SetChS2Sel(0.);     // Center of channel selection window
+  tofAnaTestbeam->SetDChS2Sel(100.);  // Width  of channel selection window
+  tofAnaTestbeam->SetTShift(0.);      // Shift DTD4 to 0
+  tofAnaTestbeam->SetSel2TOff(0.);    // Shift Sel2 time peak to 0
+  tofAnaTestbeam->SetTOffD4(13.);     // Shift DTD4 to physical value
+
+  tofAnaTestbeam->SetBeamRefSmType(iRSelTyp);  // common reaction reference
+  tofAnaTestbeam->SetBeamRefSmId(iRSelSm);
+  tofAnaTestbeam->SetBeamRefRpc(iRSelRpc);
+
+  if (iSel2 > -1) {
+    tofClust->SetSel2Id(iSel2);
+    tofClust->SetSel2Sm(iSel2Sm);
+    tofClust->SetSel2Rpc(iSel2Rpc);
+
+    tofAnaTestbeam->SetMrpcSel2(
+      iSel2);  // initialization of second selector Mrpc Type
+    tofAnaTestbeam->SetMrpcSel2Sm(
+      iSel2Sm);  // initialization of second selector Mrpc SmId
+    tofAnaTestbeam->SetMrpcSel2Rpc(
+      iSel2Rpc);  // initialization of second selector Mrpc RpcId
+  }
+
+  tofClust->SetDutId(iDut);
+  tofClust->SetDutSm(iDutSm);
+  tofClust->SetDutRpc(iDutRpc);
+
+  tofClust->SetSelId(iRef);
+  tofClust->SetSelSm(iRefSm);
+  tofClust->SetSelRpc(iRefRpc);
+
+  tofAnaTestbeam->SetDut(iDut);            // Device under test
+  tofAnaTestbeam->SetDutSm(iDutSm);        // Device under test
+  tofAnaTestbeam->SetDutRpc(iDutRpc);      // Device under test
+  tofAnaTestbeam->SetMrpcRef(iRef);        // Reference RPC
+  tofAnaTestbeam->SetMrpcRefSm(iRefSm);    // Reference RPC
+  tofAnaTestbeam->SetMrpcRefRpc(iRefRpc);  // Reference RPC
+
+  tofAnaTestbeam->SetChi2Lim(10.);  // initialization of Chi2 selection limit
+  cout << "Run with iRSel = " << iRSel << ", iSel2 = " << iSel2in << endl;
+
+  if (0) switch (iSet) {
+      case 0:   // upper part of setup: P2 - P5
+      case 3:   // upper part of setup: P2 - P5
+      case 34:  // upper part of setup: P2 - P5
+      case 400300:
+        switch (iRSel) {
+          case 4:
+            tofAnaTestbeam->SetTShift(0.);    // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(16.);   // Shift DTD4 to physical value
+            tofAnaTestbeam->SetSel2TOff(0.);  // Shift Sel2 time peak to 0
+            break;
+
+          case 5:
+            tofAnaTestbeam->SetTShift(-3.);   // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(16.);   // Shift DTD4 to physical value
+            tofAnaTestbeam->SetSel2TOff(0.);  // Shift Sel2 time peak to 0
+            break;
+
+          case 9:
+            tofAnaTestbeam->SetChi2Lim(
+              100.);  // initialization of Chi2 selection limit
+            tofAnaTestbeam->SetMulDMax(
+              3);  // Max Multiplicity in BeamRef // Diamond
+            tofAnaTestbeam->SetTShift(0.1);    // Shift DTD4 to 0
+            tofAnaTestbeam->SetTOffD4(16.);    // Shift DTD4 to physical value
+            tofAnaTestbeam->SetSel2TOff(0.5);  // Shift Sel2 time peak to 0
+            break;
+
+          default:;
+        }
+
+      default:
+        cout << "<E> detector setup " << iSet << " unknown, stop!" << endl;
+        return;
+        ;
+    }  // end of different subsets
+
+  run->AddTask(tofAnaTestbeam);
+  // =========================================================================
+  /*
+   CbmTofOnlineDisplay* display = new CbmTofOnlineDisplay();
+   display->SetUpdateInterval(1000);
+   run->AddTask(display);   
+   */
+  // -----  Parameter database   --------------------------------------------
+
+  FairRuntimeDb* rtdb       = run->GetRuntimeDb();
+  Bool_t kParameterMerged   = kTRUE;
+  FairParRootFileIo* parIo2 = new FairParRootFileIo(kParameterMerged);
+  parIo2->open(ParFile.Data(), "UPDATE");
+  parIo2->print();
+  rtdb->setFirstInput(parIo2);
+
+  FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
+  parIo1->open(parFileList, "in");
+  parIo1->print();
+  rtdb->setSecondInput(parIo1);
+  rtdb->print();
+  rtdb->printParamContexts();
+
+  //  FairParRootFileIo* parInput1 = new FairParRootFileIo();
+  //  parInput1->open(ParFile.Data());
+  //  rtdb->setFirstInput(parInput1);
+
+  FairEventManager* fMan = new FairEventManager();
+
+  CbmEvDisTracks* Tracks = new CbmEvDisTracks(
+    "Tof Tracks",
+    1,
+    kFALSE,
+    kTRUE);  //name, verbosity, RnrChildren points, RnrChildren track
+             //  CbmEvDisTracks *Tracks =  new CbmEvDisTracks("Tof Tracks",1);
+  fMan->AddTask(Tracks);
+  CbmPixelHitSetDraw* TofUHits =
+    new CbmPixelHitSetDraw("TofUHit", kRed, kOpenCross);
+  fMan->AddTask(TofUHits);
+  CbmPointSetArrayDraw* TofHits = new CbmPointSetArrayDraw(
+    "TofHit",
+    1,
+    1,
+    4,
+    kTRUE);  //name, colorMode, markerMode, verbosity, RnrChildren
+  //  CbmPixelHitSetDraw *TofHits = new CbmPixelHitSetDraw ("TofHit", kRed, kOpenCircle, 4);// kFullSquare);
+  fMan->AddTask(TofHits);
+
+
+  TGeoVolume* top = gGeoManager->GetTopVolume();
+  gGeoManager->SetVisOption(1);
+  gGeoManager->SetVisLevel(5);
+  TObjArray* allvolumes = gGeoManager->GetListOfVolumes();
+  //cout<<"GeoVolumes  "  << gGeoManager->GetListOfVolumes()->GetEntries()<<endl;
+  for (Int_t i = 0; i < allvolumes->GetEntries(); i++) {
+    TGeoVolume* vol = (TGeoVolume*) allvolumes->At(i);
+    TString name    = vol->GetName();
+    //    cout << " GeoVolume "<<i<<" Name: "<< name << endl;
+    vol->SetTransparency(90);
+    /* switch (char *) not allowed any more in root 6 :(
+    switch(name.Data()) {
+    case "counter":
+      vol->SetTransparency(95);
+      break;
+
+    case "tof_glass":
+    case "Gap":
+    case "Cell":
+      vol->SetTransparency(99);
+      break;
+
+    case "pcb":
+      vol->SetTransparency(30);
+      break;
+
+    default:
+      vol->SetTransparency(96);
+    }
+    */
+  }
+  //  gGeoManager->SetVisLevel(3);
+  //  top->SetTransparency(80);
+  //  top->Draw("ogl");
+
+  //  fMan->Init(1,4,10000);
+  fMan->Init(1, 5);
+
+  cout << "customize TEveManager gEve " << gEve << endl;
+  gEve->GetDefaultGLViewer()->SetClearColor(kYellow - 10);
+  TGLViewer* v       = gEve->GetDefaultGLViewer();
+  TGLAnnotation* ann = new TGLAnnotation(v, cFileId, 0.01, 0.98);
+  ann->SetTextSize(0.03);  // % of window diagonal
+  ann->SetTextColor(4);
+
+  //  gEve->TEveProjectionAxes()->SetDrawOrigin(kTRUE);
+
+  {  // from readCurrentCamera(const char* fname)
+    TGLCamera& c      = gEve->GetDefaultGLViewer()->CurrentCamera();
+    const char* fname = "Cam.sav";
+    TFile* f          = TFile::Open(fname, "READ");
+    if (!f) return;
+    if (f->GetKey(c.ClassName())) {
+      f->GetKey(c.ClassName())->Read(&c);
+      c.IncTimeStamp();
+      gEve->GetDefaultGLViewer()->RequestDraw();
+    }
+  }
+
+  // -----   Intialise and run   --------------------------------------------
+  // run->Init();
+  //  cout << "Starting run" << endl;
+  //  run->Run(0, nEvents);
+  // ------------------------------------------------------------------------
+  // default display
+  /*
+  TString Display_Status = "pl_over_Mat04D4best.C";
+  TString Display_Funct = "pl_over_Mat04D4best()";  
+  gROOT->LoadMacro(Display_Status);
+
+  gROOT->LoadMacro("fit_ybox.h");
+  gROOT->LoadMacro("pl_all_CluMul.C");
+  gROOT->LoadMacro("pl_all_CluRate.C");
+  gROOT->LoadMacro("pl_over_cluSel.C");
+  gROOT->LoadMacro("pl_over_clu.C");
+  gROOT->LoadMacro("pl_all_dTSel.C");
+  gROOT->LoadMacro("pl_over_MatD4sel.C");
+  gROOT->LoadMacro("save_hst.C");
+
+  switch(iSet){
+    default:
+  case 0:
+  case 3:
+  case 49:
+  case 79:
+  case 34:
+  case 94:
+  case 37:
+  case 97:
+  case 39:
+  case 99:
+  case 93:
+  case 300400:
+  case 400300:
+  case 910900:
+  case 300900:
+  case 400900:
+  case 901900:
+  case 921920:
+  case 300921:
+  case 920921:
+  case 920300:
+  case 921300:
+
+    gInterpreter->ProcessLine("pl_over_clu(6)");
+    gInterpreter->ProcessLine("pl_over_clu(6,0,1)");
+    gInterpreter->ProcessLine("pl_over_clu(9,0,0)");
+    gInterpreter->ProcessLine("pl_over_clu(9,0,1)");
+    gInterpreter->ProcessLine("pl_over_clu(9,1,0)");
+    gInterpreter->ProcessLine("pl_over_clu(9,1,1)");
+    gInterpreter->ProcessLine("pl_over_clu(9,2,0)");
+    gInterpreter->ProcessLine("pl_over_clu(9,2,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,6,0,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,6,0,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,0,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,0,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,1,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,1,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,2,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(0,9,2,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,6,0,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,6,0,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,0,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,0,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,1,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,1,1)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,2,0)");
+    gInterpreter->ProcessLine("pl_over_cluSel(1,9,2,1)");
+    gInterpreter->ProcessLine("pl_all_CluMul()");
+    gInterpreter->ProcessLine("pl_all_CluRate()");
+    gInterpreter->ProcessLine("pl_all_dTSel()");
+    TString FSave=Form("save_hst(\"cosdev-status%d_%d_Cal_%s.hst.root\")",iCalSet,iSel2in,cCalId.Data());
+    gInterpreter->ProcessLine(FSave.Data());
+    //gInterpreter->ProcessLine("pl_over_MatD4sel()");
+    break;
+    ;
+  }
+  gInterpreter->ProcessLine(Display_Funct);
+  */
+}
diff --git a/macro/beamtime/mcbm2018/dis_trks.C b/macro/beamtime/mcbm2018/dis_trks.C
new file mode 100644
index 0000000000000000000000000000000000000000..c3e06005cce5616c2f80bbc90195e6857e81c8b1
--- /dev/null
+++ b/macro/beamtime/mcbm2018/dis_trks.C
@@ -0,0 +1,652 @@
+void dis_trks(Int_t nEvents        = 10,
+              Int_t iSel           = 1,
+              Int_t iGenCor        = 1,
+              TString cFileId      = "48.50.7.1",
+              TString cSet         = "000010020",
+              Int_t iSel2          = 20,
+              Int_t iTrackingSetup = 2,
+              Double_t dScalFac    = 1.,
+              Double_t dChi2Lim2   = 500.,
+              Double_t dDeadtime   = 50,
+              TString cCalId       = "",
+              Int_t iAnaCor        = 1,
+              Bool_t bUseSigCalib  = kFALSE,
+              Int_t iCalSet        = 30040500,
+              Int_t iCalOpt        = 1,
+              Int_t iMc            = 0) {
+
+  Int_t iVerbose = 1;
+  if (cCalId == "") cCalId = cFileId;
+  TString FId = cFileId;
+  TString cRun(FId(0, 3));
+  Int_t iRun = cRun.Atoi();
+  // Specify log level (INFO, DEBUG, DEBUG1, ...)
+  //TString logLevel = "FATAL";
+  //TString logLevel = "ERROR";
+  //TString logLevel = "INFO";
+  TString logLevel = "DEBUG";
+  //TString logLevel = "DEBUG1";
+  //TString logLevel = "DEBUG2";
+  //TString logLevel = "DEBUG3";
+  TString workDir = gSystem->Getenv("VMCWORKDIR");
+  //TString workDir          = "../../..";
+  //TString paramDir       = workDir  + "/macro/beamtime/mcbm2018";
+  TString paramDir      = "./";
+  TString ParFile       = paramDir + "/data/" + cFileId.Data() + ".params.root";
+  TString InputFile     = paramDir + "/data/" + cFileId.Data() + ".root";
+  TString InputDigiFile = paramDir + "/data/digidev_" + cFileId.Data()
+                          + Form("_%s_%02.0f_Cal", cSet.Data(), dDeadtime)
+                          + cCalId + ".out.root";
+  if (iMc == 1) {
+    InputFile     = paramDir + "/data/" + cFileId.Data() + ".raw.root";
+    InputDigiFile = paramDir + "/data/" + cFileId.Data() + ".rec.root";
+    iRun          = 700;
+  }
+  TString OutputFile = paramDir + "/data/hits_" + cFileId.Data()
+                       + Form("_%s_%06d_%03d", cSet.Data(), iSel, iSel2)
+                       + ".out.root";
+  TString cHstFile =
+    paramDir
+    + Form(
+        "/hst/%s_%03.0f_%s_%06d_%03d_%03.1f_%03.1f_trk%03d_Cal%s_Dis.hst.root",
+        cFileId.Data(),
+        dDeadtime,
+        cSet.Data(),
+        iSel,
+        iSel2,
+        dScalFac,
+        dChi2Lim2,
+        iTrackingSetup,
+        cCalId.Data());
+  TString cTrkFile = Form("%s_tofFindTracks.hst.root", cCalId.Data());
+  TString cAnaFile = Form("%s_TrkAnaTestBeam.hst.root", cCalId.Data());
+
+  cout << " InputDigiFile = " << InputDigiFile << endl;
+
+  TString shcmd = "rm -v " + ParFile;
+  gSystem->Exec(shcmd.Data());
+
+  TList* parFileList = new TList();
+
+  Int_t iGeo = 0;  //iMc;
+  if (iGeo == 0) {
+    TString TofGeo = "v18m_mcbm";
+    cout << "Geometry version " << TofGeo << endl;
+
+    // -----   Load the geometry setup   -------------------------------------
+    /*
+  const char* setupName = "mcbm_beam_2020_03";
+  TString setupFile = workDir + "/geometry/setup/setup_" + setupName + ".C";
+  TString setupFunct = "setup_";
+  setupFunct = setupFunct + setupName + "()";
+  std::cout << "-I- mcbm_reco: Loading macro " << setupFile << std::endl;
+  gROOT->LoadMacro(setupFile);
+  gROOT->ProcessLine(setupFunct);
+  CbmSetup* setup = CbmSetup::Instance();
+*/
+
+    TString geoDir  = gSystem->Getenv("VMCWORKDIR");
+    TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+    TFile* fgeo     = new TFile(geoFile);
+    TGeoManager* geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+    if (NULL == geoMan) {
+      cout << "<E> FAIRGeom not found in geoFile" << endl;
+      return;
+    }
+
+    //TObjString *tofDigiFile = new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par"); // TOF digi file
+    //parFileList->Add(tofDigiFile);
+
+    TObjString* tofDigiBdfFile =
+      new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digibdf.par");
+    parFileList->Add(tofDigiBdfFile);
+
+    // -----   Reconstruction run   -------------------------------------------
+    FairRunAna* run = new FairRunAna();
+    cout << "InputFile:     " << InputFile.Data() << endl;
+    cout << "InputDigiFile: " << InputDigiFile.Data() << endl;
+
+    //run->SetInputFile(InputFile.Data());
+    //run->AddFriend(InputDigiFile.Data());
+    run->SetInputFile(InputDigiFile.Data());
+    //run->AddFriend(InputFile.Data());
+    run->SetOutputFile(OutputFile);
+
+    FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+    //  FairLogger::GetLogger()->SetLogVerbosityLevel("MEDIUM");
+    FairLogger::GetLogger()->SetLogVerbosityLevel("VERYHIGH");
+
+    // -----   Local selection variables  -------------------------------------------
+
+    Int_t iRef    = iSel % 1000;
+    Int_t iDut    = (iSel - iRef) / 1000;
+    Int_t iDutRpc = iDut % 10;
+    iDut          = (iDut - iDutRpc) / 10;
+    Int_t iDutSm  = iDut % 10;
+    iDut          = (iDut - iDutSm) / 10;
+    Int_t iRefRpc = iRef % 10;
+    iRef          = (iRef - iRefRpc) / 10;
+    Int_t iRefSm  = iRef % 10;
+    iRef          = (iRef - iRefSm) / 10;
+
+    Int_t iSel2in  = iSel2;
+    Int_t iSel2Rpc = iSel2 % 10;
+    iSel2          = (iSel2 - iSel2Rpc) / 10;
+    Int_t iSel2Sm  = iSel2 % 10;
+    iSel2          = (iSel2 - iSel2Sm) / 10;
+
+    Int_t calMode = 93;
+    Int_t calSel  = 1;
+    Bool_t bOut   = kFALSE;
+
+    CbmTofEventClusterizer* tofClust =
+      new CbmTofEventClusterizer("TOF Event Clusterizer", iVerbose, bOut);
+    Int_t calSelRead = calSel;
+    if (calSel < 0) calSelRead = 0;
+    TString cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                          cFileId.Data(),
+                          iCalSet,
+                          calMode,
+                          calSelRead);
+    if (cCalId != "XXX")
+      cFname = Form("%s_set%09d_%02d_%01dtofClust.hst.root",
+                    cCalId.Data(),
+                    iCalSet,
+                    calMode,
+                    calSelRead);
+    tofClust->SetCalParFileName(cFname);
+    TString cOutFname =
+      Form("tofClust_%s_set%09d.hst.root", cFileId.Data(), iCalSet);
+    tofClust->SetOutHstFileName(cOutFname);
+
+    // =========================================================================
+    // ===                       Tracking                                    ===
+    // =========================================================================
+
+    CbmTofTrackFinder* tofTrackFinder = new CbmTofTrackFinderNN();
+    tofTrackFinder->SetMaxTofTimeDifference(0.4);  // in ns/cm
+    Int_t TrackerPar = 0;
+    switch (TrackerPar) {
+      case 0:                           // for full mTof setup
+        tofTrackFinder->SetTxLIM(0.3);  // max slope dx/dz
+        tofTrackFinder->SetTyLIM(0.3);  // max dev from mean slope dy/dz
+        tofTrackFinder->SetTxMean(0.);  // mean slope dy/dz
+        tofTrackFinder->SetTyMean(0.);  // mean slope dy/dz
+        break;
+      case 1:                             // for double stack test counters
+        tofTrackFinder->SetTxMean(0.21);  // mean slope dy/dz
+        tofTrackFinder->SetTyMean(0.18);  // mean slope dy/dz
+        tofTrackFinder->SetTxLIM(0.15);   // max slope dx/dz
+        tofTrackFinder->SetTyLIM(0.18);   // max dev from mean slope dy/dz
+        break;
+    }
+
+    CbmTofTrackFitter* tofTrackFitter = new CbmTofTrackFitterKF(0, 211);
+    TFitter* MyFit                    = new TFitter(1);  // initialize Minuit
+    tofTrackFinder->SetFitter(tofTrackFitter);
+    CbmTofFindTracks* tofFindTracks = new CbmTofFindTracks("TOF Track Finder");
+    tofFindTracks->UseFinder(tofTrackFinder);
+    tofFindTracks->UseFitter(tofTrackFitter);
+    tofFindTracks->SetCalOpt(
+      iCalOpt);  // 1 - update offsets, 2 - update walk, 0 - bypass
+    tofFindTracks->SetCorMode(
+      iGenCor);  // valid options: 0,1,2,3,4,5,6, 10 - 19
+    tofFindTracks->SetTtTarg(
+      0.057);  // target value for inverse velocity, > 0.033 ns/cm!
+    //tofFindTracks->SetTtTarg(0.035);            // target value for inverse velocity, > 0.033 ns/cm!
+    tofFindTracks->SetCalParFileName(
+      cTrkFile);  // Tracker parameter value file name
+    tofFindTracks->SetBeamCounter(5, 0, 0);  // default beam counter
+    tofFindTracks->SetR0Lim(100.);
+
+    tofFindTracks->SetStationMaxHMul(
+      30);  // Max Hit Multiplicity in any used station
+
+    tofFindTracks->SetT0MAX(dScalFac);  // in ns
+    tofFindTracks->SetSIGT(0.08);       // default in ns
+    tofFindTracks->SetSIGX(0.3);        // default in cm
+    tofFindTracks->SetSIGY(0.6);        // default in cm
+    tofFindTracks->SetSIGZ(0.05);       // default in cm
+    tofFindTracks->SetUseSigCalib(
+      bUseSigCalib);  // ignore resolutions in CalPar file
+    tofTrackFinder->SetSIGLIM(dChi2Lim2
+                              * 2.);  // matching window in multiples of chi2
+    tofTrackFinder->SetChiMaxAccept(dChi2Lim2);  // max tracklet chi2
+
+
+    Int_t iMinNofHits   = -1;
+    Int_t iNStations    = 0;
+    Int_t iNReqStations = 3;
+
+    switch (iTrackingSetup) {
+      case 0:  // bypass mode
+        iMinNofHits = -1;
+        iNStations  = 1;
+        tofFindTracks->SetStation(0, 5, 0, 0);  // Diamond
+        break;
+
+      case 1:  // for calibration mode of full setup
+      {
+        Double_t dTsig = dScalFac * 0.03;
+        tofFindTracks->SetSIGT(dTsig);  // allow for variable deviations in ns
+      }
+        iMinNofHits   = 3;
+        iNStations    = 28;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 2, 2);
+        tofFindTracks->SetStation(2, 0, 1, 2);
+        tofFindTracks->SetStation(3, 0, 0, 2);
+        tofFindTracks->SetStation(4, 0, 2, 1);
+        tofFindTracks->SetStation(5, 0, 1, 1);
+        tofFindTracks->SetStation(6, 0, 0, 1);
+        tofFindTracks->SetStation(7, 0, 2, 3);
+        tofFindTracks->SetStation(8, 0, 1, 3);
+        tofFindTracks->SetStation(9, 0, 0, 3);
+        tofFindTracks->SetStation(10, 0, 2, 0);
+        tofFindTracks->SetStation(11, 0, 1, 0);
+        tofFindTracks->SetStation(12, 0, 0, 0);
+        tofFindTracks->SetStation(13, 0, 2, 4);
+        tofFindTracks->SetStation(14, 0, 1, 4);
+        tofFindTracks->SetStation(15, 0, 0, 4);
+        tofFindTracks->SetStation(16, 0, 4, 0);
+        tofFindTracks->SetStation(17, 0, 3, 0);
+        tofFindTracks->SetStation(18, 0, 4, 1);
+        tofFindTracks->SetStation(19, 0, 3, 1);
+        tofFindTracks->SetStation(20, 0, 4, 2);
+        tofFindTracks->SetStation(21, 0, 3, 2);
+        tofFindTracks->SetStation(22, 0, 4, 3);
+        tofFindTracks->SetStation(23, 0, 3, 3);
+        tofFindTracks->SetStation(24, 0, 4, 4);
+        tofFindTracks->SetStation(25, 0, 3, 4);
+        tofFindTracks->SetStation(26, 9, 0, 0);
+        tofFindTracks->SetStation(27, 9, 0, 1);
+        //tofFindTracks->SetStation(28, 6, 0, 0);
+        //tofFindTracks->SetStation(29, 6, 0, 1);
+        break;
+
+      case 2:
+        iMinNofHits   = 3;
+        iNStations    = 14;
+        iNReqStations = 5;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 1);
+        tofFindTracks->SetStation(2, 0, 3, 1);
+        tofFindTracks->SetStation(3, 0, 4, 0);
+        tofFindTracks->SetStation(4, 0, 3, 0);
+        tofFindTracks->SetStation(5, 0, 4, 2);
+        tofFindTracks->SetStation(6, 0, 3, 2);
+        tofFindTracks->SetStation(7, 0, 4, 3);
+        tofFindTracks->SetStation(8, 0, 3, 3);
+        tofFindTracks->SetStation(9, 0, 4, 4);
+        tofFindTracks->SetStation(10, 0, 3, 4);
+        tofFindTracks->SetStation(11, 9, 0, 0);
+        tofFindTracks->SetStation(12, 9, 0, 1);
+        tofFindTracks->SetStation(13, 7, 0, 0);
+        break;
+
+      case 3:
+        iMinNofHits   = 3;
+        iNStations    = 16;
+        iNReqStations = 4;
+        
+		tofFindTracks->SetStation(0, 0, 2, 2);
+		tofFindTracks->SetStation(1, 0, 1, 2);
+		tofFindTracks->SetStation(2, 0, 0, 2);
+
+		tofFindTracks->SetStation(3, 0, 2, 1);
+		tofFindTracks->SetStation(4, 0, 1, 1);
+		tofFindTracks->SetStation(5, 0, 0, 1);
+
+		tofFindTracks->SetStation(6, 0, 2, 3);
+		tofFindTracks->SetStation(7, 0, 1, 3);
+		tofFindTracks->SetStation(8, 0, 0, 3);
+
+		tofFindTracks->SetStation(9,  0, 2, 0);
+		tofFindTracks->SetStation(10, 0, 1, 0);
+		tofFindTracks->SetStation(11, 0, 0, 0);
+
+		tofFindTracks->SetStation(12, 0, 2, 4);
+		tofFindTracks->SetStation(13, 0, 1, 4);
+		tofFindTracks->SetStation(14, 0, 0, 4);
+
+		tofFindTracks->SetStation(15, 5, 0, 0);
+        /*
+     tofFindTracks->SetStation(16, 0, 3, 2);         
+     tofFindTracks->SetStation(17, 0, 4, 2);  
+     tofFindTracks->SetStation(18, 0, 3, 1);         
+     tofFindTracks->SetStation(19, 0, 4, 1);
+     tofFindTracks->SetStation(20, 0, 3, 3);         
+     tofFindTracks->SetStation(21, 0, 4, 3);
+     tofFindTracks->SetStation(22, 0, 3, 0);         
+     tofFindTracks->SetStation(23, 0, 4, 0);
+     tofFindTracks->SetStation(24, 0, 3, 4);         
+     tofFindTracks->SetStation(25, 0, 4, 4); 
+     */
+        break;
+
+      case 4:
+        iMinNofHits   = 3;
+        iNStations    = 4;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 2, 2);
+        tofFindTracks->SetStation(2, 0, 0, 2);
+        tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+        break;
+
+      case 5:  // for calibration of 2-stack and add-on counters (STAR2, BUC)
+        iMinNofHits   = 3;
+        iNStations    = 5;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 1);
+        tofFindTracks->SetStation(2, 0, 3, 1);
+        tofFindTracks->SetStation(3, 9, 0, 0);
+        tofFindTracks->SetStation(4, 9, 0, 1);
+        break;
+
+      case 6:
+        iMinNofHits   = 3;
+        iNStations    = 4;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 1);
+        tofFindTracks->SetStation(2, 0, 3, 1);
+        tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+        //    tofFindTracks->SetStation(3, 9, 0, 0);
+        //    tofFindTracks->SetStation(3, 9, 0, 1);
+        //    tofFindTracks->SetStation(3, 7, 0, 0);
+        break;
+
+      case 7:  // for calibration of 2-stack and add-on counters (BUC)
+        iMinNofHits   = 4;
+        iNStations    = 5;
+        iNReqStations = 5;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 3);
+        tofFindTracks->SetStation(2, 0, 3, 3);
+        tofFindTracks->SetStation(3, 6, 0, 0);
+        tofFindTracks->SetStation(4, 6, 0, 1);
+        break;
+
+      case 8:  // evaluation of add-on counters (BUC)
+        iMinNofHits   = 3;
+        iNStations    = 4;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 4, 3);
+        tofFindTracks->SetStation(2, 0, 3, 3);
+        tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+        break;
+
+      case 10:
+        iMinNofHits   = 3;
+        iNStations    = 4;
+        iNReqStations = 4;
+        tofFindTracks->SetStation(0, 5, 0, 0);
+        tofFindTracks->SetStation(1, 0, 1, 2);
+        tofFindTracks->SetStation(2, 0, 0, 2);
+        tofFindTracks->SetStation(3, 0, 2, 2);
+
+      default:
+        cout << "Tracking setup " << iTrackingSetup << " not implemented "
+             << endl;
+        return;
+        ;
+    }
+    tofFindTracks->SetMinNofHits(iMinNofHits);
+    tofFindTracks->SetNStations(iNStations);
+    tofFindTracks->SetNReqStations(iNReqStations);
+    tofFindTracks->PrintSetup();
+    run->AddTask(tofFindTracks);
+
+    // =========================================================================
+    // ===                       Analysis                                    ===
+    // =========================================================================
+    CbmTofAnaTestbeam* tofAnaTestbeam =
+      new CbmTofAnaTestbeam("TOF TestBeam Analysis", iVerbose);
+    tofAnaTestbeam->SetCorMode(iAnaCor);  // 1 - DTD4, 2 - X4, 3 - Y4, 4 - Texp
+    tofAnaTestbeam->SetHitDistMin(30.);   // initialization
+    tofAnaTestbeam->SetEnableMatchPosScaling(kTRUE);
+    tofAnaTestbeam->SetSpillDuration(4.);
+    //CbmTofAnaTestbeam defaults
+    tofAnaTestbeam->SetDXMean(0.);
+    tofAnaTestbeam->SetDYMean(0.);
+    tofAnaTestbeam->SetDTMean(0.);  // in ns
+    tofAnaTestbeam->SetDXWidth(0.5);
+    tofAnaTestbeam->SetDYWidth(1.0);
+    tofAnaTestbeam->SetDTWidth(0.1);  // in ns
+    tofAnaTestbeam->SetCalParFileName(cAnaFile);
+    Double_t dScalFacA = 0.9;  // dScalFac is used for tracking
+    tofAnaTestbeam->SetPosY4Sel(
+      0.5 * dScalFacA);  // Y Position selection in fraction of strip length
+    tofAnaTestbeam->SetDTDia(0.);    // Time difference to additional diamond
+    tofAnaTestbeam->SetMul0Max(20);  // Max Multiplicity in dut
+    tofAnaTestbeam->SetMul4Max(30);  // Max Multiplicity in Ref - RPC
+    tofAnaTestbeam->SetMulDMax(3);   // Max Multiplicity in Diamond / BeamRef
+    tofAnaTestbeam->SetTOffD4(14.);  // initialization
+    tofAnaTestbeam->SetDTD4MAX(
+      6.);  // initialization of Max time difference Ref - BRef
+
+    //tofAnaTestbeam->SetTShift(-28000.);// initialization
+    tofAnaTestbeam->SetPosYS2Sel(
+      0.55);  // Y Position selection in fraction of strip length
+    tofAnaTestbeam->SetChS2Sel(0.);     // Center of channel selection window
+    tofAnaTestbeam->SetDChS2Sel(100.);  // Width  of channel selection window
+    tofAnaTestbeam->SetSel2TOff(0.);    // Shift Sel2 time peak to 0
+    tofAnaTestbeam->SetChi2Lim(5.);  // initialization of Chi2 selection limit
+    tofAnaTestbeam->SetChi2Lim2(
+      3.);  // initialization of Chi2 selection limit for Mref-Sel2 pair
+    tofAnaTestbeam->SetDutDX(
+      15.);  // limit inspection of tracklets to selected region
+    tofAnaTestbeam->SetDutDY(
+      15.);  // limit inspection of tracklets to selected region
+    tofAnaTestbeam->SetSIGLIM(3.);  // max matching chi2
+    tofAnaTestbeam->SetSIGT(0.08);  // in ns
+    tofAnaTestbeam->SetSIGX(0.3);   // in cm
+    tofAnaTestbeam->SetSIGY(0.6);   // in cm
+
+    Int_t iRSel    = 500;
+    Int_t iRSelTyp = 5;
+    Int_t iRSelSm  = 0;
+    Int_t iRSelRpc = 0;
+    Int_t iRSelin  = iRSel;
+
+    tofAnaTestbeam->SetBeamRefSmType(iRSelTyp);  // common reaction reference
+    tofAnaTestbeam->SetBeamRefSmId(iRSelSm);
+    tofAnaTestbeam->SetBeamRefRpc(iRSelRpc);
+
+    if (iSel2 >= 0) {
+      tofAnaTestbeam->SetMrpcSel2(
+        iSel2);  // initialization of second selector Mrpc Type
+      tofAnaTestbeam->SetMrpcSel2Sm(
+        iSel2Sm);  // initialization of second selector Mrpc SmId
+      tofAnaTestbeam->SetMrpcSel2Rpc(
+        iSel2Rpc);  // initialization of second selector Mrpc RpcId
+    }
+
+    tofAnaTestbeam->SetDut(iDut);            // Device under test
+    tofAnaTestbeam->SetDutSm(iDutSm);        // Device under test
+    tofAnaTestbeam->SetDutRpc(iDutRpc);      // Device under test
+    tofAnaTestbeam->SetMrpcRef(iRef);        // Reference RPC
+    tofAnaTestbeam->SetMrpcRefSm(iRefSm);    // Reference RPC
+    tofAnaTestbeam->SetMrpcRefRpc(iRefRpc);  // Reference RPC
+
+    cout << "dispatch iSel = " << iSel << ", iSel2in = " << iSel2in
+         << ", iRSelin = " << iRSelin << ", iRSel = " << iRSel << endl;
+    if (1) {
+      switch (iSel) {
+
+        case 10:
+          switch (iRSelin) {
+            case 500:
+              tofAnaTestbeam->SetTShift(2.5);  // Shift DTD4 to 0
+              tofAnaTestbeam->SetTOffD4(18.);  // Shift DTD4 to physical value
+              switch (iSel2in) {
+                case 20:
+                  tofAnaTestbeam->SetSel2TOff(0.);  // Shift Sel2 time peak to 0
+                  break;
+                default:;
+              }
+              break;
+            default:;
+          }
+          break;
+
+        case 700040:
+        case 900040:
+        case 901040:
+          switch (iRSelin) {
+            case 500:
+              tofAnaTestbeam->SetTShift(0.3);  // Shift DTD4 to 0
+              tofAnaTestbeam->SetTOffD4(18.);  // Shift DTD4 to physical value
+
+              switch (iSel2in) {
+                case 30:
+                  tofAnaTestbeam->SetSel2TOff(
+                    -0.3);  // Shift Sel2 time peak to 0
+                  break;
+                case 31:
+                  tofAnaTestbeam->SetSel2TOff(
+                    -0.41);  // Shift Sel2 time peak to 0
+                  break;
+
+                default:;
+              }
+              break;
+            default:;
+          }
+          break;
+
+        case 700041:
+        case 900041:
+        case 901041:
+          switch (iRSelin) {
+            case 500:
+              tofAnaTestbeam->SetTShift(5.3);  // Shift DTD4 to 0
+              tofAnaTestbeam->SetTOffD4(18.);  // Shift DTD4 to physical value
+
+              switch (iSel2in) {
+                case 30:
+                  tofAnaTestbeam->SetSel2TOff(
+                    -0.3);  // Shift Sel2 time peak to 0
+                  break;
+                case 31:
+                  tofAnaTestbeam->SetSel2TOff(
+                    -0.41);  // Shift Sel2 time peak to 0
+                  break;
+
+                default:;
+              }
+              break;
+            default:;
+          }
+          break;
+
+        case 600043:
+        case 601043:
+          switch (iRSelin) {
+            case 500:
+              tofAnaTestbeam->SetTShift(5.3);  // Shift DTD4 to 0
+              tofAnaTestbeam->SetTOffD4(11.);  // Shift DTD4 to physical value
+
+              switch (iSel2in) {
+                case 33:
+                  tofAnaTestbeam->SetSel2TOff(
+                    -0.55);  // Shift Sel2 time peak to 0
+                  break;
+
+                default:;
+              }
+              break;
+            default:;
+          }
+          break;
+
+        default:
+          cout
+            << "Better to define analysis setup! Running with default offset "
+               "parameter... "
+            << endl;
+          // return;
+      }  // end of different subsets
+
+      cout << " Initialize TSHIFT to " << tofAnaTestbeam->GetTShift() << endl;
+      //run->AddTask(tofAnaTestbeam);
+    }
+
+    // -----  Parameter database   --------------------------------------------
+    FairRuntimeDb* rtdb       = run->GetRuntimeDb();
+    Bool_t kParameterMerged   = kTRUE;
+    FairParRootFileIo* parIo2 = new FairParRootFileIo(kParameterMerged);
+    parIo2->open(ParFile.Data(), "UPDATE");
+    parIo2->print();
+    rtdb->setFirstInput(parIo2);
+
+    FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
+    parIo1->open(parFileList, "in");
+    parIo1->print();
+    rtdb->setSecondInput(parIo1);
+    rtdb->print();
+    rtdb->printParamContexts();
+
+    //  FairParRootFileIo* parInput1 = new FairParRootFileIo();
+    //  parInput1->open(ParFile.Data());
+    //  rtdb->setFirstInput(parInput1);
+
+    FairEventManager* fMan = new FairEventManager();
+
+    CbmEvDisTracks* Tracks = new CbmEvDisTracks(
+      "Tof Tracks",
+      1,
+      kFALSE,
+      kTRUE);  //name, verbosity, RnrChildren points, RnrChildren track
+               //  CbmEvDisTracks *Tracks =  new CbmEvDisTracks("Tof Tracks",1);
+    fMan->AddTask(Tracks);
+    CbmPixelHitSetDraw* TofUHits =
+      new CbmPixelHitSetDraw("TofUHit", kRed, kOpenCross);
+    fMan->AddTask(TofUHits);
+    CbmPointSetArrayDraw* TofHits = new CbmPointSetArrayDraw(
+      "TofHit",
+      1,
+      1,
+      1,
+      kTRUE);  //name, colorMode, markerMode, verbosity, RnrChildren
+    //  CbmPixelHitSetDraw *TofHits = new CbmPixelHitSetDraw ("TofHit", kRed, kOpenCircle, 4);// kFullSquare);
+    fMan->AddTask(TofHits);
+
+    TGeoVolume* top = gGeoManager->GetTopVolume();
+    gGeoManager->SetVisOption(1);
+    gGeoManager->SetVisLevel(5);
+    TObjArray* allvolumes = gGeoManager->GetListOfVolumes();
+    //cout<<"GeoVolumes  "  << gGeoManager->GetListOfVolumes()->GetEntries()<<endl;
+    for (Int_t i = 0; i < allvolumes->GetEntries(); i++) {
+      TGeoVolume* vol = (TGeoVolume*) allvolumes->At(i);
+      //TString name = vol->GetName();
+      //cout << " GeoVolume "<<i<<" Name: "<< name << endl;
+      vol->SetLineColor(kRed);
+      vol->SetTransparency(80);
+    }
+    fMan->Init(1, 5);
+
+    cout << "customize TEveManager gEve " << gEve << endl;
+    gEve->GetDefaultGLViewer()->SetClearColor(kYellow - 10);
+    TGLViewer* v       = gEve->GetDefaultGLViewer();
+    TGLAnnotation* ann = new TGLAnnotation(v, cFileId, 0.01, 0.98);
+    ann->SetTextSize(0.03);  // % of window diagonal
+    ann->SetTextColor(4);
+
+    {  // from readCurrentCamera(const char* fname)
+      TGLCamera& c      = gEve->GetDefaultGLViewer()->CurrentCamera();
+      const char* fname = "Cam.sav";
+      TFile* f          = TFile::Open(fname, "READ");
+      if (!f) return;
+      if (f->GetKey(c.ClassName())) {
+        f->GetKey(c.ClassName())->Read(&c);
+        c.IncTimeStamp();
+        gEve->GetDefaultGLViewer()->RequestDraw();
+      }
+    }
+  }
+}
diff --git a/macro/beamtime/mcbm2018/init_cal_all.sh b/macro/beamtime/mcbm2018/init_cal_all.sh
new file mode 100755
index 0000000000000000000000000000000000000000..eab6b0a83ee5a4478ac42a60c499f441ef96e8b7
--- /dev/null
+++ b/macro/beamtime/mcbm2018/init_cal_all.sh
@@ -0,0 +1,146 @@
+#!/bin/bash
+# shell script to initialize clusterizer calibrations
+#SBATCH -J calall
+#SBATCH -D /lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018
+#SBATCH --time=8:00:00
+##SBATCH --time=6-00:00:00
+#SBATCH --mem=2000
+##SBATCH --partition=long
+cRun=$1
+
+echo 'Initialize clusterizer calibration for run '$cRun
+
+iCalSet=$2
+((iTmp  = $iCalSet ))
+((iBRef = $iTmp % 1000))
+((iTmp  = $iTmp - $iBRef))
+((iSet  = $iTmp / 1000))
+((iMRef = $iTmp % 1000000))
+((iMRef = $iMRef / 1000))
+((iTmp  = $iTmp - $iMRef))
+((iDut  = $iTmp / 1000000))
+echo Calib setup is ${iCalSet}, iSet=$iSet, iDut=$iDut, iMRef=$iMRef, iBRef=$iBRef
+cCalSet=$iCalSet
+if (( iCalSet<100000000 )); then 
+cCalSet="0"$iCalSet
+fi
+if (( iCalSet<10000000 )); then 
+cCalSet="00"$iCalSet
+fi
+if (( iCalSet<1000000 )); then 
+cCalSet="000"$iCalSet
+fi
+if (( iCalSet<100000 )); then 
+cCalSet="0000"$iCalSet
+fi
+echo cCalSet = $cCalSet
+#iSet=0
+#lastOpt=''
+nEvi0=150000 # start value
+nEvi1=50000  # increment 
+
+if [ -e /lustre ]; then
+source /lustre/cbm/users/nh/CBM/cbmroot/trunk/build/config.sh 
+wdir=/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018
+outdir=/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/${cRun}
+else 
+wdir=`pwd`
+outdir=${wdir}/${cRun}
+fi
+mkdir ${outdir}
+
+cd  ${wdir}
+mkdir ${cRun}
+cp rootlogon.C ${cRun}
+cp .rootrc ${cRun}
+cd ${cRun}
+
+# Global variables, for for-loops
+iRestart=0
+#iRestart=13
+iStep=0
+iStepLast=0
+iCalSel0=0
+iCalSel1=1
+#iCalSel0=-3 #0
+#iCalSel1=-4 #1
+# ************************** Starting while Loop ***************************** #
+(( nEvi = nEvi0 + 10*nEvi1 ))
+optList=""
+optList=`echo " $nEvi,93,1,$iMRef,0 "`$optList 
+icalmod=3
+for icallev in 8 8 7 7 6 5 4 4 3 3 1
+do
+    (( nEvi = nEvi0 + (icallev-1)*nEvi1 ))
+    optList=`echo " $nEvi,$icallev$icalmod,$iCalSel0,$iDut,0 "`$optList
+    if [ $iMRef -ne 14 ]; then 
+	optList=`echo " $nEvi,$icallev$icalmod,$iCalSel1,$iMRef,0 "`$optList
+    else 
+	for iMod in 40  10 
+	do
+	    if [ $iMod -ne $iDut ]; then
+		optList=`echo " $nEvi,$icallev$icalmod,$iCalSel1,$iMod,0 "`$optList
+	    fi
+	done
+    fi
+    if [ $icallev -lt 7 ]; then
+      optList=`echo " $nEvi,$icallev$icalmod,$iCalSel0,$iBRef,50 "`$optList 
+      optList=`echo " $nEvi,$icallev$icalmod,$iCalSel1,$iBRef,50 "`$optList
+    else
+#      optList=`echo " $nEvi,$icallev$icalmod,-2,2,0 "`$optList
+      echo skip additional step
+    fi
+done
+ optList=`echo " $nEvi,0,0,$iBRef,50 "`$optList      # start Init1
+ echo optList:  $optList
+
+for inOpt in $optList
+do  
+    echo step ${iStep} with option $inOpt
+    ((iStepLast = ${iStep}))
+    ((iStep += 1))
+
+    mkdir Init${iStep}
+    cp rootlogon.C Init${iStep}
+    cp .rootrc Init${iStep}
+    cd Init${iStep}
+
+    if [[ ${lastOpt:+1} ]] ; then
+	# echo last round was done with $lastOpt, extract 2. and 3. word
+	i1=`expr index $inOpt , `
+	i2=($i1+3)
+	#echo `expr index $inOpt , ` = $i1
+	cMode=${inOpt:$i1:2}
+	cSel=${inOpt:$i2:1}
+	echo next iteration: cMode=$cMode, cSel=$cSel 
+	if [[ ${cSel} = "-" ]];then 
+	    cSel=${inOpt:$i2:2}
+	    echo cSel=$cSel 
+	    cSel="0"
+	fi
+	#copy calibration file 
+	if (($iStep > $iRestart)) ; then
+	    cp -v ../Init${iStepLast}/tofClust_${cRun}_set${cCalSet}.hst.root ${cRun}_set${cCalSet}_${cMode}_${cSel}tofClust.hst.root
+	fi
+    fi 
+
+    lastOpt=$inOpt
+    # generate new calibration file
+    if (($iStep > $iRestart)) ; then 
+	root -b -q '../../ana_digi_cal_all.C('$inOpt',"'${cRun}'",'${iCalSet}',0,'${iBRef}') '
+
+	cp -v tofClust_${cRun}_set${cCalSet}.hst.root ../${cRun}_set${cCalSet}_${cMode}_${cSel}tofClust.hst.root
+	cp *pdf ../
+	#./screenshot.sh
+	cd .. 
+	rm ../${cRun}_set${cCalSet}_${cMode}_${cSel}tofClust.hst.root
+	ln -s ./${cRun}/${cRun}_set${cCalSet}_${cMode}_${cSel}tofClust.hst.root ../${cRun}_set${cCalSet}_${cMode}_${cSel}tofClust.hst.root
+	echo Init step $iStep with mode ${cMode}, option $inOpt  finished
+    else 
+	cd ..
+	echo Init step $iStep with mode ${cMode}, option $inOpt  skipped
+    fi   
+done
+
+cd  ${wdir}
+mv -v slurm-${SLURM_JOB_ID}.out ${outdir}/InitCalib_${cRun}_${cCalSet}.out
diff --git a/macro/beamtime/mcbm2018/pl_all_Track2D.C b/macro/beamtime/mcbm2018/pl_all_Track2D.C
new file mode 100644
index 0000000000000000000000000000000000000000..3d9fc5dd2eb5ddf4c7ff75fd9851ac57260f0fd0
--- /dev/null
+++ b/macro/beamtime/mcbm2018/pl_all_Track2D.C
@@ -0,0 +1,122 @@
+void pl_all_Track2D(Int_t iOpt = 1, Int_t iNSt = 4) {
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2);
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas* can = new TCanvas("can", "can", 48, 56, 900, 900);
+  can->Divide(5, 6, 0.01, 0.01);
+  //  can->Divide(2,2,0,0);
+  Float_t lsize = 0.07;
+
+  gPad->SetFillColor(0);
+  gStyle->SetPalette(1);
+  gStyle->SetLabelSize(lsize);
+
+  //gStyle->SetOptStat(kTRUE);
+  //gROOT->cd();
+  //gROOT->SetDirLevel(2);
+
+  TH2* h;
+  TH2* h2;
+  const Int_t iType[6]   = {0, 9, 6, 5, 7, 8};
+  const Int_t iSmNum[6]  = {5, 1, 1, 1, 1, 1};
+  const Int_t iRpcNum[6] = {5, 2, 2, 1, 1, 8};
+  TString cOpt;
+
+  switch (iOpt) {
+    case 0: cOpt = "Size"; break;
+    case 1: cOpt = "Pos"; break;
+    case 2: cOpt = "TOff"; break;
+    case 3: cOpt = "Tot"; break;
+    case 4: cOpt = "Walk"; break;
+    case 5: cOpt = "Walk"; break;
+    case 6: cOpt = "Mul"; break;
+    case 7: cOpt = "Trms"; break;
+    case 8: cOpt = "DelPos"; break;
+    case 9: cOpt = "DelTOff"; break;
+    case 10: cOpt = "DelMatPos"; break;
+    case 11: cOpt = "DelMatTOff"; break;
+    default:;
+  }
+
+  Int_t iDet       = 0;
+  Double_t dAvMean = 0.;
+  Double_t dAvRMS  = 0.;
+  Int_t iCanv      = 0;
+
+  for (Int_t iSt = 0; iSt < iNSt; iSt++) {
+    // cout << "plot station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+    for (Int_t iSm = 0; iSm < iSmNum[iSt]; iSm++) {
+      //cout << "plot module at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+      for (Int_t iRp = 0; iRp < iRpcNum[iSt]; iRp++) {
+        //cout << "plot rpc at station "<<iSt<<" with "<< iSmNum[iSt] <<" modules of "<<iRpcNum[iSt]<<" Rpcs each"<<endl;
+        can->cd(iCanv + 1);
+        iCanv++;
+        gROOT->cd();
+        TString hname = "";
+        Int_t iCol    = 1;
+        switch (iOpt) {
+          case 4:
+            for (Int_t iSide = 0; iSide < 2; iSide++)
+              for (Int_t iCh = 0; iCh < 32; iCh++) {
+                hname = Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_%s",
+                             iType[iSt],
+                             iSm,
+                             iRp,
+                             iCh,
+                             iSide,
+                             cOpt.Data());
+                h     = (TH2*) gROOT->FindObjectAny(hname);
+                if (h != NULL) {
+                  TProfile* hProf =
+                    h->ProfileX(Form("%s_pfx%d%d", hname.Data(), iCh, iSide));
+                  hProf->SetLineColor(iCol);
+                  hProf->SetLineStyle(1);
+                  hProf->SetMarkerColor(iCol);
+                  hProf->SetMarkerStyle(24 + iSide);
+                  iCol++;
+                  if (iCh == 0) iCol = 1;
+                  if (iCh == 0 && iSide == 0) {
+                    hProf->SetMaximum(0.4);
+                    hProf->SetMinimum(-0.4);
+                    hProf->GetXaxis()->SetRangeUser(0., 10.);
+                    hProf->Draw("LP");
+                  } else {
+                    hProf->Draw("LPsame");
+                  }
+                }
+              }
+            break;
+          default:
+            hname = Form("cal_SmT%01d_sm%03d_rpc%03d_%s",
+                         iType[iSt],
+                         iSm,
+                         iRp,
+                         cOpt.Data());
+            h     = (TH2*) gROOT->FindObjectAny(hname);
+            if (h != NULL) {
+              if (iOpt == 2 || iOpt == 2) { gPad->SetLogz(); }
+              h->Draw("colz");
+              h->ProfileX()->Draw("same");
+              iDet++;
+              dAvMean += h->ProfileX()->GetMean(2);
+              dAvRMS += h->ProfileX()->GetRMS(2);
+              cout << "TrackQA " << cOpt.Data() << " for TSR " << iType[iSt]
+                   << iSm << iRp << ": Off " << h->ProfileX()->GetMean(2)
+                   << ", RMS " << h->ProfileX()->GetRMS(2) << endl;
+            }
+        }
+      }
+    }
+  }
+  dAvMean /= (Double_t) iDet;
+  dAvRMS /= (Double_t) iDet;
+  cout << "TrackQA " << cOpt.Data() << ": AvOff " << dAvMean << ", AvRMS "
+       << dAvRMS << endl;
+  dAvMean = TMath::Abs(dAvMean);
+  gROOT->ProcessLine(
+    Form(".! echo %d > %sAvOff.res", (Int_t)(dAvMean * 1.E4), cOpt.Data()));
+  gROOT->ProcessLine(
+    Form(".! echo %d > %sAvRMS.res", (Int_t)(dAvRMS * 1.E4), cOpt.Data()));
+
+  can->SaveAs(Form("pl_all_Track_%s.pdf", cOpt.Data()));
+}
diff --git a/macro/beamtime/mcbm2018/trk_cal_digi.sh b/macro/beamtime/mcbm2018/trk_cal_digi.sh
new file mode 100755
index 0000000000000000000000000000000000000000..6917f26ef4146e8df85af246c914c136b45ea911
--- /dev/null
+++ b/macro/beamtime/mcbm2018/trk_cal_digi.sh
@@ -0,0 +1,251 @@
+#!/bin/bash
+# shell script to apply clusterizer calibrations
+#SBATCH -J trk_cal_digi
+#SBATCH -D /lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018
+##SBATCH --time=8:00:00
+#SBATCH --time=5-24:00:00
+#SBATCH --mem=4000
+#SBATCH --partition=long
+
+trk_cal_digi() {
+cRun=$1
+
+iCalSet=$2
+((iTmp  = $iCalSet ))
+((iBRef = $iTmp % 1000))
+((iTmp  = $iTmp - $iBRef))
+((iSet  = $iTmp / 1000))
+((iRef  = $iTmp % 1000000))
+((iRef  = $iRef / 1000))
+((iTmp  = $iTmp - $iRef))
+((iDut  = $iTmp / 1000000))
+
+iSel2=$3
+cSel2=$iSel2
+if (( iSel2<100 )); then 
+cSel2="0"$iSel2
+fi
+if (( iSel2<10 )); then 
+cSel2="00"$iSel2
+fi
+
+cCalSet=$iCalSet
+if (( iCalSet<100000000 )); then 
+cCalSet="0"$iCalSet
+fi
+if (( iCalSet<10000000 )); then 
+cCalSet="00"$iCalSet
+fi
+if (( iCalSet<1000000 )); then 
+cCalSet="000"$iCalSet
+fi
+if (( iCalSet<100000 )); then 
+cCalSet="0000"$iCalSet
+fi
+echo cCalSet = $cCalSet
+
+Deadtime=$4
+if [[ ${Deadtime} = "" ]]; then
+Deadtime=50.
+fi
+
+CalIdMode=$5
+if [[ ${CalIdMode} = "" ]]; then
+ echo use native calibration file 
+ CalIdMode=${cRun}
+ CalFile=${cRun}_set${cCalSet}_93_1tofClust.hst.root
+else 
+ CalFile=${CalIdMode}_set${cCalSet}_93_1tofClust.hst.root
+ RunFile=${cRun}_set${cCalSet}_93_1tofClust.hst.root
+# rm ${RunFile}
+# ln -s ${CalFile} ${RunFile} 
+ echo use calibrations from  ${CalFile}
+fi
+
+iCalOpt=$6
+if [[ ${iCalOpt} = "" ]]; then
+  iCalOpt=1	
+fi
+
+iTraSetup=$7
+if [[ $iTraSetup = "" ]]; then 
+  iTraSetup=1
+fi
+
+CalIdSet=$8
+if [[ ${CalIdSet} = "" ]]; then
+    echo use native calibration file
+    CalIdSet=$cCalSet
+else
+    CalFile=${CalIdMode}_set${CalIdSet}_93_1tofClust.hst.root    
+fi
+
+
+echo trk_cal_digi for $cRun with iDut=$iDut, iRef=$iRef, iSet=$iCalSet, iSel2=$iSel2, iBRef=$iBRef, Deadtime=$Deadtime, CalFile=$CalFile
+
+if [[ $iShLev = "" ]]; then 
+  iShLev=0
+  nEvt=200000
+  dDTres=100000
+  dDTRMSres=100000
+fi 
+
+echo execute trk_cal_digi at shell level $iShLev
+
+if [ -e /lustre/cbm ]; then
+source /lustre/cbm/users/nh/CBM/cbmroot/trunk/build/config.sh 
+wdir=/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018
+outdir=/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/${cRun}
+else 
+wdir=`pwd`
+outdir=${wdir}/${cRun}
+fi
+
+iter=0;
+
+cd $wdir
+mkdir $cRun
+cd    $cRun 
+cp    ../.rootrc .
+cp    ../rootlogon.C .
+
+echo Execute in `pwd` at $iShLev: ./trk_cal_digi.sh $1 $2 $3 $4 $5 $6 $7 $8
+
+# get initial digi calibration 
+#cp -v  ./I*/${CalFile}  .
+ 
+# get latest tracker offsets
+# cp -v ../${cRun}_tofFindTracks.hst.root .
+
+rm -v TCalib.res
+nEvtMax=0
+(( nEvtMax = nEvt*10 ))
+
+#frange1 limits DT spectrum range 
+fRange1=2.
+# frange2 limits chi2
+fRange2=4.0
+TRange2Limit=3. 
+
+#iSel=911921
+iSel=10020
+iGenCor=3
+cCalSet2=${cCalSet}_$cSel2
+
+case $iCalOpt in
+  1)   # TOff
+  ;;
+  2)   # Walk
+  (( nEvt *= 10 ))
+  (( nEvtMax *= 10 ))
+  ;;
+esac
+
+while [[ $dDTres > 0 ]]; do
+  nEvt=`echo "scale=0;$nEvt * 1./1." | bc`
+  #nEvt=`echo "scale=0;$nEvt * 1.1/1." | bc`
+  if [ $nEvt -gt $nEvtMax ]; then
+    nEvt=$nEvtMax
+  fi
+
+  fRange2=`echo "$fRange2 * 0.9" | bc`
+  compare_TRange2=`echo "$fRange2 < $TRange2Limit" | bc`
+  if  [[ $compare_TRange2 > 0 ]]; then
+   fRange2=$TRange2Limit
+  fi
+  
+  iCalAct=$iCalOpt
+  iIter=0
+  echo Enter while loop with CalAct $iCalAct in dir `pwd`
+  while [[ $iCalAct -gt 0 ]]; do  
+    cd $wdir/$cRun
+    echo Current loop with CalAct $iCalAct and CalOpt $iCalOpt
+    if [[ $iCalOpt = 1 ]] || [[ $iCalAct -gt 1 ]]; then 
+      root -b -q '../ana_digi_cal.C('$nEvt',93,1,'$iRef',1,"'$cRun'",'$iCalSet',1,'$iSel2','$Deadtime',"'$CalIdMode'") '
+      # update calibration parameter file, will only be active in next iteration 
+      if [[ $iIter = -10 ]] && [[ $iCalOpt = 1 ]]; then  # exploratory option when iIter set to 0 
+        echo Update Calibration file from ana_digi_cal
+        cp -v tofClust_${cRun}_set${cCalSet}.hst.root ../${cRun}_set${cCalSet}_93_1tofClust.hst.root
+        echo 20000 > TOffAvOff.res
+        echo 20000 > TOffAvRMS.res
+      else
+        root -b -q '../ana_trks.C('$nEvt','$iSel','$iGenCor',"'$cRun'","'$cCalSet2'",'$iSel2','$iTraSetup','$fRange1','$fRange2','$Deadtime',"'$CalIdMode'",1,1,'$iCalSet','$iCalAct')'
+      #root -l 'ana_trksi.C(-1,10,1,"385.50.5.0","000014500_020",20,1,1.90,7.60,50,"385.50.5.0",1,1)'
+  
+        cp -v New_${CalFile} ${CalFile}  
+      fi
+    else 
+      cd $wdir
+      # store current status 
+      dLDTres=$dDTres
+      dLDTRMSres=$dDTRMSres
+      iLCalOpt=$iCalOpt
+      echo Store limits $dLDTres, $dLDTRMSres
+      (( iShLev += 1 ))
+      echo exec in `pwd` at level $iShLev: trk_cal_digi $1 $2 $3 $4 $5 1 $7
+      trk_cal_digi $1 $2 $3 $4 $5 1 $7
+      (( iShLev -= 1 ))
+      # restore old status
+      dDTres=$dLDTres
+      dDTRMSres=$dLDTRMSres
+      iCalOpt=$iLCalOpt
+      echo exec1done, resume old CalOpt $iCalOpt with status $dDTres, $dDTRMSres
+    fi
+    (( iCalAct -= 1 ))
+    (( iIter   += 1 ))
+    echo Continue while loop with CalAct $iCalAct and CalOpt $iCalOpt
+  done
+  
+  cd $wdir/$cRun
+  Tres=`cat TOffAvOff.res`
+  TRMSres=`cat TOffAvRMS.res`
+
+  if [[ $Tres = 0 ]]; then
+    Tres=1
+  fi
+  dTdif=`echo "$dDTres - $Tres" | bc`
+  compare_result=`echo "$Tres < $dDTres" | bc`
+
+  dTRMSdif=`echo "$dDTRMSres - $TRMSres" | bc`
+  compare_RMS=`echo "$TRMSres < $dDTRMSres" | bc`
+
+  echo at iter=$iter got TOff = $Tres, compare to $dDTres, dTdif = $dTdif, result = $compare_result, TRMS = $TRMSres, old $dDTRMSres, dif = $dTRMSdif, result = $compare_RMS 
+
+  ((compare_result += $compare_RMS))
+  echo CMPR result_summary: $compare_result 
+
+#  if [ $iter = 1 ]; then 
+#    exit 0  # for debugging 
+#  fi
+
+  if [[ $compare_result > 0 ]]; then
+    if [[ $Tres = 0 ]]; then
+      Tres=1
+    fi
+    dDTres=$Tres
+    dDTRMSres=$TRMSres
+    echo Store new res values $dDTres, $dDTRMSres
+    (( dDTRMSres -= 1 ))  # next attempt should be at least 1ps better for continuation
+    cp -v New_${CalFile} ${CalFile}  
+    cp -v New_${CalFile} ${CalFile}_$iter  
+  else
+    dDTres=0
+  fi
+  (( iter += 1 ))
+done
+
+cd $wdir/$cRun
+# generate full statistics digi file 
+if [[ $iShLev = 0 ]]; then 
+  root -b -q '../ana_digi_cal.C(-1,93,1,'$iRef',1,"'$cRun'",'$iCalSet',1,'$iSel2','$Deadtime',"'$CalIdMode'") '
+fi
+
+cd $wdir
+
+if [[ $iShLev = 0 ]]; then 
+  mv -v slurm-${SLURM_JOB_ID}.out ${outdir}/TrkCalDigi_${cRun}_${iCalSet}_${iSel2}_${CalIdMode}.out
+fi
+
+} #end of function body
+
+trk_cal_digi $1 $2 $3 $4 $5 $6 $7 $8
diff --git a/macro/beamtime/pl_Track2D.C b/macro/beamtime/pl_Track2D.C
new file mode 100644
index 0000000000000000000000000000000000000000..73cde6d7ac7855cf58ec0f19d1c6cfce3f83890f
--- /dev/null
+++ b/macro/beamtime/pl_Track2D.C
@@ -0,0 +1,129 @@
+void pl_Track2D(Int_t iOpt = 1, Int_t iCounterId=22, Int_t iStrip=-1, Double_t TotMax=10.) {
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas* can = new TCanvas("can", "can", 48, 56, 900, 900);
+  //can->Divide(4, 4, 0.01, 0.01);
+  //  can->Divide(2,2,0,0);
+  if(iStrip > -1) can->Divide(1,2,0.01,0.01);
+  
+  Float_t lsize = 0.07;
+
+  gPad->SetFillColor(0);
+  gStyle->SetPalette(1);
+  gStyle->SetLabelSize(lsize);
+
+  //gStyle->SetOptStat(kTRUE);
+  //gROOT->cd();
+  //gROOT->SetDirLevel(2);
+
+  TH2* h;
+  TH2* h2;
+
+  TString cOpt;
+
+  switch (iOpt) {
+    case 0: cOpt = "Size"; break;
+    case 1: cOpt = "Pos"; break;
+    case 2: cOpt = "TOff"; break;
+    case 3: cOpt = "Tot"; break;
+    case 4: 
+    case 14:
+		    cOpt = "Walk"; break;
+    case 5: cOpt = "Walk"; break;
+    case 6: cOpt = "Mul"; break;
+    case 7: cOpt = "Trms"; break;
+    case 8: cOpt = "DelPos"; break;
+    case 9: cOpt = "DelTOff"; break;
+    case 10: cOpt = "DelMatPos"; break;
+    case 11: cOpt = "DelMatTOff"; break;
+    default:;
+  }
+
+  Int_t iDet       = 0;
+  Double_t dAvMean = 0.;
+  Double_t dAvRMS  = 0.;
+  Int_t iCanv      = 0;
+
+  Int_t iRp=iCounterId%10;
+  Int_t iSm=((iCounterId-iRp)/10) % 10;
+  Int_t iType=(iCounterId-iSm*10-iRp)/100;
+
+
+  gROOT->cd();
+  TString hname = "";
+  Int_t iCol    = 1;
+  switch (iOpt) {
+    case 4:
+    {
+    Int_t iDraw=0;
+    for (Int_t iSide = 0; iSide < 2; iSide++) {
+	  iCol=1;
+      for (Int_t iCh = 0; iCh < 32; iCh++) {
+		if(iStrip > -1) if (iCh != iStrip) continue;
+		cout << "CS "<<iCh<<iSide<<endl;
+        hname = Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_%s",
+                     iType,iSm,iRp,iCh,iSide,cOpt.Data());
+        h     = (TH2*) gROOT->FindObjectAny(hname);
+        if (h != NULL) {
+          TProfile* hProf =
+          h->ProfileX(Form("%s_pfx%d%d", hname.Data(), iCh, iSide));
+          hProf->SetLineColor(iCol);
+          hProf->SetLineStyle(1);
+          hProf->SetMarkerColor(iCol);
+          hProf->SetMarkerStyle(24 + iSide);
+          iCol++;
+          if (iCh == 0) iCol = 1;
+          if ( iDraw == 0 ) {
+            hProf->SetMaximum(0.4);
+            hProf->SetMinimum(-0.4);
+            hProf->GetXaxis()->SetRangeUser(0., TotMax);
+            hProf->Draw("LP");
+          } else {
+            hProf->Draw("LPsame");
+          }
+          iDraw++;
+        }
+      }
+    }
+    }
+    break;
+    
+    case 14:
+    {
+      Int_t iCanv=1;
+      for (Int_t iSide = 0; iSide < 2; iSide++)
+        for (Int_t iCh = 0; iCh < 32; iCh++) {
+		  if(iStrip > -1) if (iCh != iStrip) continue;
+          hname = Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_%s",
+                     iType,iSm,iRp,iCh,iSide,cOpt.Data());
+          h     = (TH2*) gROOT->FindObjectAny(hname);
+          if (h != NULL) {
+            can->cd(iCanv++);
+            h->Draw("colz");
+            gPad->SetLogz();
+            h->ProfileX()->Draw("same");
+            TString hName2 = Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_Walk_px",
+                                  iType,iSm,iRp,iCh,iSide);
+            TH1* hCor = (TH1*) gROOT->FindObjectAny(hName2);
+            if ( NULL != hCor ) {
+			  hCor->SetLineColor(6);
+			  hCor->Draw("same");
+			}                      
+          }
+        }
+    }
+    break;
+        
+    default:
+      hname = Form("cal_SmT%01d_sm%03d_rpc%03d_%s",
+                   iType,iSm,iRp,cOpt.Data());
+      h     = (TH2*) gROOT->FindObjectAny(hname);
+      if (h != NULL) {
+        if (iOpt == 2 || iOpt == 2) { gPad->SetLogz(); }
+          h->Draw("colz");
+          h->ProfileX()->Draw("same");
+      }
+  }
+  
+  can->SaveAs(Form("pl_Track_%s_%d.pdf", cOpt.Data(),iCounterId));
+}
diff --git a/macro/beamtime/pl_trk_Walk.C b/macro/beamtime/pl_trk_Walk.C
new file mode 100644
index 0000000000000000000000000000000000000000..4597de1dc2242ba36d7a436c67db949b6ad9bdb5
--- /dev/null
+++ b/macro/beamtime/pl_trk_Walk.C
@@ -0,0 +1,133 @@
+void pl_trk_Walk(Int_t iId=900, Int_t iOpt=1, Double_t dMax=0.)
+{
+  //  TCanvas *can = new TCanvas("can22","can22");
+  //  can->Divide(2,2); 
+  //  TCanvas *can = new TCanvas("can","can",48,55,700,900);
+  TCanvas *can = new TCanvas("can","can",48,56,1400,900);
+  //  can->Divide(2,2,0,0); 
+  Float_t lsize=0.07;
+
+ gPad->SetFillColor(0);
+ gStyle->SetPalette(1);
+ gStyle->SetLabelSize(lsize);
+
+ //gStyle->SetOptStat(kTRUE);
+ //gROOT->cd();
+ //gROOT->SetDirLevel(2);
+
+ TH2 *h;
+ TH2 *h2;
+ Int_t iRp=iId%10;
+ Int_t iSm=((iId-iRp)/10 )%10;
+ Int_t iTy=(iId-iSm*10-iRp)/100;
+
+ TString cOpt="Walk";
+
+ Int_t iDet=0;
+ Int_t iCanv=0;
+
+ gROOT->cd();
+ TString hname="";
+ Int_t iCol=1;
+ switch(iOpt) {
+	case 0:;
+	  break;
+    case 1:
+      {
+	  can->Divide(4,8,0.01,0.01); 
+	  for (Int_t iCh=0; iCh<32; iCh++) {
+        can->cd(iCanv+1); iCanv++;
+	    for (Int_t iSide=0; iSide<2; iSide++) {
+		  hname=Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_%s",iTy,iSm,iRp,iCh,iSide,cOpt.Data());
+          h=(TH2 *)gROOT->FindObjectAny(hname);
+          if (h!=NULL) {
+            TProfile *hProf=h->ProfileX(Form("%s_pfx%d%d",hname.Data(),iCh,iSide));
+		    iCol=iSide+2;
+            hProf->SetLineColor(iCol);
+            hProf->SetMarkerColor(iCol);
+            hProf->SetMarkerStyle(20+iSide);
+            if(iSide==0) {
+			  if(dMax>0.) {
+			    hProf->SetMaximum(dMax);
+			    hProf->SetMinimum(-dMax);
+			  }
+			  hProf->Draw();
+			} else {
+			  hProf->Draw("same");
+			}
+		  }
+              
+		}
+	  }
+		break;
+		
+    case 20:
+      {
+	  can->Divide(4,8,0.01,0.01); 
+	  for (Int_t iCh=0; iCh<32; iCh++) {
+        can->cd(iCanv+1); iCanv++;
+	    Int_t iSide=0; 
+		hname=Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_%s",iTy,iSm,iRp,iCh,iSide,cOpt.Data());
+        h=(TH2 *)gROOT->FindObjectAny(hname);
+        if (h!=NULL) {
+			if(dMax>0.) {
+//			  h->GetYAxis()->SetMaximum(dMax);
+//			  h->GetYAxis()->SetMinimum(-dMax);
+			}
+			h->Draw("colz");
+		}
+	  }			  
+	  }
+	  break;
+	  
+    case 21:
+      {
+	  can->Divide(4,8,0.01,0.01); 
+	  for (Int_t iCh=0; iCh<32; iCh++) {
+        can->cd(iCanv+1); iCanv++;
+	    Int_t iSide=1;  
+		hname=Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_%s",iTy,iSm,iRp,iCh,iSide,cOpt.Data());
+        h=(TH2 *)gROOT->FindObjectAny(hname);
+        if (h!=NULL) {
+			h->Draw("colz");
+		}
+	  }			  
+	  }	  
+	  break;
+	  
+    case 200:
+      {
+	  can->Divide(4,8,0.01,0.01); 
+	  for (Int_t iCh=0; iCh<32; iCh++) {
+        can->cd(iCanv+1); iCanv++;
+	    Int_t iSide=0; 
+		hname=Form("cl_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_%s",iTy,iSm,iRp,iCh,iSide,cOpt.Data());
+        h=(TH2 *)gROOT->FindObjectAny(hname);
+        if (h!=NULL) {
+			h->Draw("colz");
+		}
+	  }			  
+	  }
+	  break;	  
+	  
+    case 201:
+      {
+	  can->Divide(4,8,0.01,0.01); 
+	  for (Int_t iCh=0; iCh<32; iCh++) {
+        can->cd(iCanv+1); iCanv++;
+	    Int_t iSide=1; 
+		hname=Form("cl_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_%s",iTy,iSm,iRp,iCh,iSide,cOpt.Data());
+        h=(TH2 *)gROOT->FindObjectAny(hname);
+        if (h!=NULL) {
+			h->Draw("colz");
+		}
+	  }			  
+	  }
+	  break;	
+	  
+ 		default:
+			;
+	}
+ }
+ can->SaveAs(Form("pl_trk_Walk_%d.pdf",iId));
+} 
diff --git a/macro/mcbm/geometry/tof/Create_TOF_Geometry_v20f_mcbm.C b/macro/mcbm/geometry/tof/Create_TOF_Geometry_v20f_mcbm.C
new file mode 100644
index 0000000000000000000000000000000000000000..e9336a8d66a081307d336e349e134dc62ebf1a65
--- /dev/null
+++ b/macro/mcbm/geometry/tof/Create_TOF_Geometry_v20f_mcbm.C
@@ -0,0 +1,1451 @@
+///
+/// \file Create_TOF_Geometry_v20f_mcbm.C
+/// \brief Generates TOF geometry in Root format.
+///
+
+// Changelog
+//
+// 2020-04-14 - v20b - NH - swapped double stack layer 2 with STAR2 moodule, buc kept as dummy
+// 2020-04-01 - v20a - NH - move mTOF +20 cm in x direction for the Mar 2020 run
+// 2019-11-28 - v19b - DE - move mTOF +12 cm in x direction for the Nov 2019 run
+// 2019-07-31 - v19a - DE - this TOF March 2019 geometry is also known as v18m
+// 2017-11-03 - v18i - DE - shift mTOF to z=298 cm for acceptance matching with mSTS
+// 2017-10-06 - v18h - DE - put v18f into vertical position to fit into the mCBM cave
+// 2017-07-15 - v18g - DE - swap the z-position of TOF modules: 2 in the front, 3 in the back
+// 2017-07-14 - v18f - DE - reduce vertical gap between TOF modules to fix the gap between modules 1-2 and 4-5
+// 2017-05-17 - v18e - DE - rotate electronics away from beam, shift 16 cm away from beam along x-axis
+// 2017-05-17 - v18d - DE - change geometry name to v18d
+
+// in root all sizes are given in cm
+
+#include "TFile.h"
+#include "TGeoCompositeShape.h"
+#include "TGeoManager.h"
+#include "TGeoMaterial.h"
+#include "TGeoMatrix.h"
+#include "TGeoMedium.h"
+#include "TGeoPgon.h"
+#include "TGeoVolume.h"
+#include "TList.h"
+#include "TMath.h"
+#include "TROOT.h"
+#include "TString.h"
+#include "TSystem.h"
+
+#include <iostream>
+
+// Name of geometry version and output file
+const TString geoVersion      = "tof_v20f_mcbm";  // do not change
+const TString geoVersionStand = geoVersion + "Stand";
+//
+const TString fileTag      = "tof_v20f";
+const TString FileNameSim  = fileTag + "_mcbm.geo.root";
+const TString FileNameGeo  = fileTag + "_mcbm_geo.root";
+const TString FileNameInfo = fileTag + "_mcbm.geo.info";
+
+// TOF_Z_Front corresponds to front cover of outer super module towers
+const Float_t TOF_Z_Front_Stand = 247.2;  // = z=298 mCBM@SIS18
+const Float_t TOF_X_Front_Stand = 0.;  // = z=298 mCBM@SIS18
+const Float_t TOF_Z_Front       = 0.;      // = z=298 mCBM@SIS18
+//const Float_t TOF_Z_Front =  130;  // = z=225 mCBM@SIS18
+//const Float_t TOF_Z_Front =  250;  // SIS 100 hadron
+//const Float_t TOF_Z_Front =  450;  // SIS 100 hadron
+//const Float_t TOF_Z_Front =  600;  // SIS 100 electron
+//const Float_t TOF_Z_Front =  650;  // SIS 100 muon
+//const Float_t TOF_Z_Front =  880;  // SIS 300 electron
+//const Float_t TOF_Z_Front = 1020;  // SIS 300 muon
+//
+//const Float_t TOF_Z_Front = 951.5;   // Wall_Z_Position = 1050 cm
+
+
+// Names of the different used materials which are used to build the modules
+// The materials are defined in the global media.geo file
+const TString KeepingVolumeMedium = "air";
+const TString BoxVolumeMedium     = "aluminium";
+const TString NoActivGasMedium    = "RPCgas_noact";
+const TString ActivGasMedium      = "RPCgas";
+const TString GlasMedium          = "RPCglass";
+const TString ElectronicsMedium   = "carbon";
+
+// Counters:
+// 0 MRPC3a
+// 1 MRPC3b
+// 2 USTC
+// 3
+// 4 Diamond
+//
+// 6 Buc 2019
+// 7 CERN 20gap
+// 8 Ceramic Pad
+const Int_t NumberOfDifferentCounterTypes = 9;
+const Float_t Glass_X[NumberOfDifferentCounterTypes] =
+  {32., 52., 32., 32., 0.2, 32., 28.8, 20., 2.4};
+const Float_t Glass_Y[NumberOfDifferentCounterTypes] =
+  {27.0, 53., 26.8, 10., 0.2, 10., 6., 20., 2.4};
+const Float_t Glass_Z[NumberOfDifferentCounterTypes] =
+  {0.1, 0.1, 0.1, 0.1, 0.01, 0.1, 0.1, 0.1, 0.1};
+
+const Float_t GasGap_X[NumberOfDifferentCounterTypes] =
+  {32., 52., 32., 32., 0.2, 32., 28.8, 20., 2.4};
+const Float_t GasGap_Y[NumberOfDifferentCounterTypes] =
+  {27.0, 53., 26.8, 10., 0.2, 10., 6., 20., 2.4};
+const Float_t GasGap_Z[NumberOfDifferentCounterTypes] =
+  {0.025, 0.025, 0.025, 0.025, 0.01, 0.02, 0.02, 0.02, 0.025};
+
+const Int_t NumberOfGaps[NumberOfDifferentCounterTypes] =
+  {8, 8, 8, 8, 1, 8, 10, 20, 4};
+//const Int_t NumberOfGaps[NumberOfDifferentCounterTypes] = {1,1,1,1}; //deb
+const Int_t NumberOfReadoutStrips[NumberOfDifferentCounterTypes] =
+  {32, 32, 32, 32, 8, 32, 32, 20, 1};
+//const Int_t NumberOfReadoutStrips[NumberOfDifferentCounterTypes] = {1,1,1,1}; //deb
+
+const Float_t SingleStackStartPosition_Z[NumberOfDifferentCounterTypes] =
+  {-0.6, -0.6, -0.6, -0.6, -0.1, -0.6, -0.6, -0.6, -1.};
+
+const Float_t Electronics_X[NumberOfDifferentCounterTypes] =
+  {34.0, 53.0, 32.0, 32., 0.3, 0.1, 28.8, 20., 0.1};
+const Float_t Electronics_Y[NumberOfDifferentCounterTypes] =
+  {5.0, 5.0, 1.0, 1., 0.1, 0.1, 1.0, 1.0, 0.1};
+const Float_t Electronics_Z[NumberOfDifferentCounterTypes] =
+  {0.3, 0.3, 0.3, 0.3, 0.1, 0.1, 0.1, 0.1, 0.1};
+
+const Int_t NofModuleTypes = 10;
+// 5 Diamond
+// 6 Buc
+// 7 CERN 20 gap
+// 8 Ceramic
+// 9 Star2
+// Aluminum box for all module types
+const Float_t Module_Size_X[NofModuleTypes] =
+  {180., 180., 180., 180., 180., 5., 40., 30., 22.5, 100.};
+const Float_t Module_Size_Y[NofModuleTypes] =
+  {49., 49., 74., 28., 18., 5., 12., 30., 11., 49.};
+const Float_t Module_Over_Y[NofModuleTypes] =
+  {11.5, 11.5, 11., 4.5, 4.5, 0., 0., 0., 0., 0.};
+const Float_t Module_Size_Z[NofModuleTypes] =
+  {11., 11., 13., 11., 11., 1., 12., 6., 6.2, 11.2};
+const Float_t Module_Thick_Alu_X_left  = 0.1;
+const Float_t Module_Thick_Alu_X_right = 1.0;
+const Float_t Module_Thick_Alu_Y       = 0.1;
+const Float_t Module_Thick_Alu_Z       = 0.1;
+
+// Distance to the center of the TOF wall [cm];
+const Float_t Wall_Z_Position = 400.;
+const Float_t MeanTheta       = 0.;
+
+//Type of Counter for module
+const Int_t CounterTypeInModule[NofModuleTypes] =
+  {0, 0, 1, 2, 3, 4, 6, 7, 8, 2};
+const Int_t NCounterInModule[NofModuleTypes] = {5, 5, 3, 5, 5, 1, 2, 1, 8, 2};
+
+// Placement of the counter inside the module
+const Float_t CounterXStartPosition[NofModuleTypes] =
+  {-60., -66.0, -56.0, -60.0, -60.0, 0.0, 0., 0., -7., 0.};
+const Float_t CounterXDistance[NofModuleTypes] =
+  {30.0, 32.0, 51.0, 30.0, 30.0, 0.0, 0., 0., 2., -1.};
+const Float_t CounterYStartPosition[NofModuleTypes] =
+  {0.0, 0.0, 0.0, 0.0, 0.0, 0., 0., -4., -1.3, 0.};
+const Float_t CounterYDistance[NofModuleTypes] =
+  {0.0, 0.0, 0.0, 0.0, 0.0, 0., 0., 8., 0., 1.};
+const Float_t CounterZDistance[NofModuleTypes] =
+  {-2.5, 0.0, 0.0, 2.5, 2.5, 0., 6., 0., 0.1, 4.};
+const Float_t CounterZStartPosition[NofModuleTypes] =
+  {0.0, 0.0, 0.0, 0.0, 0.0, 0., -3., 0., 0.0, -2.};
+const Float_t CounterRotationAngle[NofModuleTypes] =
+  {0., 8.7, 7.0, 0., 0., 0., 0., 0., 0., 0.};
+
+// Pole (support structure)
+const Int_t MaxNumberOfPoles = 20;
+Float_t Pole_ZPos[MaxNumberOfPoles];
+Float_t Pole_Col[MaxNumberOfPoles];
+Int_t NumberOfPoles = 0;
+
+const Float_t Pole_Size_X  = 20.;
+const Float_t Pole_Size_Y  = 300.;
+const Float_t Pole_Size_Z  = 10.;
+const Float_t Pole_Thick_X = 5.;
+const Float_t Pole_Thick_Y = 5.;
+const Float_t Pole_Thick_Z = 5.;
+
+// Bars (support structure)
+const Float_t Bar_Size_X = 20.;
+const Float_t Bar_Size_Y = 20.;
+Float_t Bar_Size_Z       = 100.;
+
+const Int_t MaxNumberOfBars = 20;
+Float_t Bar_ZPos[MaxNumberOfBars];
+Float_t Bar_XPos[MaxNumberOfBars];
+Int_t NumberOfBars = 0;
+
+const Float_t ChamberOverlap = 40;
+const Float_t DxColl         = 158.0;  //Module_Size_X-ChamberOverlap;
+//const Float_t Pole_Offset=Module_Size_X/2.+Pole_Size_X/2.;
+const Float_t Pole_Offset = 90.0 + Pole_Size_X / 2.;
+
+// Position for module placement
+const Float_t Inner_Module_First_Y_Position = 16.;
+const Float_t Inner_Module_Last_Y_Position  = 480.;
+const Float_t Inner_Module_X_Offset         = 0.;  // centered position in x/y
+//const Float_t Inner_Module_X_Offset=18; // shift by 16 cm in x
+const Int_t Inner_Module_NTypes                       = 3;
+const Float_t Inner_Module_Types[Inner_Module_NTypes] = {4., 3., 0.};
+//const Float_t Inner_Module_Number[Inner_Module_NTypes] = {2.,2.,6.}; //V13_3a
+const Float_t Inner_Module_Number[Inner_Module_NTypes] = {2., 2., 1.};  //V13_3a
+//const Float_t Inner_Module_Number[Inner_Module_NTypes] = {0.,0.,0.}; //debugging
+
+const Float_t InnerSide_Module_X_Offset                    = 51.;
+const Float_t InnerSide_Module_NTypes                      = 1;
+const Float_t InnerSide_Module_Types[Inner_Module_NTypes]  = {5.};
+const Float_t InnerSide_Module_Number[Inner_Module_NTypes] = {2.};  //v13_3a
+//const Float_t InnerSide_Module_Number[Inner_Module_NTypes] = {0.};  //debug
+
+const Float_t Outer_Module_First_Y_Position = 0.;
+const Float_t Outer_Module_Last_Y_Position  = 480.;
+const Float_t Outer_Module_X_Offset         = 3.;
+const Int_t Outer_Module_Col                = 4;
+const Int_t Outer_Module_NTypes             = 2;
+const Float_t Outer_Module_Types[Outer_Module_NTypes][Outer_Module_Col] =
+  {1., 1., 1., 1., 2., 2., 2., 2.};
+const Float_t Outer_Module_Number[Outer_Module_NTypes][Outer_Module_Col] =
+  {9., 9., 2., 0., 0., 0., 3., 4.};  //V13_3a
+//const Float_t Outer_Module_Number[Outer_Module_NTypes][Outer_Module_Col] = {1.,1.,0.,0.,  0.,0.,0.,0.};//debug
+
+const Float_t Star2_First_Z_Position       = TOF_Z_Front + 16.5;
+const Float_t Star2_Delta_Z_Position       = 0.;
+const Float_t Star2_First_Y_Position       = 30.35;  //
+const Float_t Star2_Delta_Y_Position       = 0.;   //
+const Float_t Star2_rotate_Z               = -90.;
+const Int_t Star2_NTypes                   = 1;
+const Float_t Star2_Types[Star2_NTypes]    = {9.};
+const Float_t Star2_Number[Star2_NTypes]   = {1.};   //debugging, V16b
+const Float_t Star2_X_Offset[Star2_NTypes] = {49.};  //{51.};
+
+const Float_t Buc_First_Z_Position     = TOF_Z_Front + 50.;
+const Float_t Buc_Delta_Z_Position     = 0.;
+const Float_t Buc_First_Y_Position     = -32.5;  //
+const Float_t Buc_Delta_Y_Position     = 0.;     //
+const Float_t Buc_rotate_Z             = 180.;
+const Int_t Buc_NTypes                 = 1;
+const Float_t Buc_Types[Buc_NTypes]    = {6.};
+const Float_t Buc_Number[Buc_NTypes]   = {1.};  //debugging, V16b
+const Float_t Buc_X_Offset[Buc_NTypes] = {53.5};
+
+const Int_t Cer_NTypes                   = 3;
+const Float_t Cer_Z_Position[Cer_NTypes] = {(float) (TOF_Z_Front + 13.2),
+                                            (float) (TOF_Z_Front + 45.),
+                                            (float) (TOF_Z_Front + 45.)};
+const Float_t Cer_X_Position[Cer_NTypes] = {0., 49.8, 49.8};
+const Float_t Cer_Y_Position[Cer_NTypes] = {-1., 5., 5.};
+const Float_t Cer_rotate_Z[Cer_NTypes]   = {0., 0., 0.};
+const Float_t Cer_Types[Cer_NTypes]      = {5., 8., 8.};
+const Float_t Cer_Number[Cer_NTypes]     = {1., 1., 1.};  //V16b
+
+const Float_t CERN_Z_Position          = TOF_Z_Front + 50;  // 20 gap
+const Float_t CERN_First_Y_Position    = 36.;
+const Float_t CERN_X_Offset            = 46.;  //65.5;
+const Float_t CERN_rotate_Z            = 90.;
+const Int_t CERN_NTypes                = 1;
+const Float_t CERN_Types[CERN_NTypes]  = {7.};  // this is the SmType!
+const Float_t CERN_Number[CERN_NTypes] = {
+  1.};  // evtl. double for split signals
+
+// some global variables
+TGeoManager* gGeoMan = NULL;           // Pointer to TGeoManager instance
+TGeoVolume* gModules[NofModuleTypes];  // Global storage for module types
+TGeoVolume* gCounter[NumberOfDifferentCounterTypes];
+TGeoVolume* gPole;
+TGeoVolume* gBar[MaxNumberOfBars];
+
+const Float_t Dia_Z_Position         = -0.2 - TOF_Z_Front_Stand;
+const Float_t Dia_First_Y_Position   = 0.;
+const Float_t Dia_X_Offset           = 0.;
+const Float_t Dia_rotate_Z           = 0.;
+const Int_t Dia_NTypes               = 1;
+const Float_t Dia_Types[Dia_NTypes]  = {5.};
+const Float_t Dia_Number[Dia_NTypes] = {1.};
+
+Float_t Last_Size_Y = 0.;
+Float_t Last_Over_Y = 0.;
+
+// Forward declarations
+void create_materials_from_media_file();
+TGeoVolume* create_counter(Int_t);
+TGeoVolume* create_new_counter(Int_t);
+TGeoVolume* create_tof_module(Int_t);
+TGeoVolume* create_new_tof_module(Int_t);
+TGeoVolume* create_tof_pole();
+TGeoVolume* create_tof_bar();
+void position_tof_poles(Int_t);
+void position_tof_bars(Int_t);
+void position_inner_tof_modules(Int_t);
+void position_side_tof_modules(Int_t);
+void position_outer_tof_modules(Int_t);
+void position_Dia(Int_t);
+void position_Star2(Int_t);
+void position_Buc(Int_t);
+void position_cer_modules(Int_t);
+void position_CERN(Int_t);
+void dump_info_file();
+
+
+void Create_TOF_Geometry_v20f_mcbm() {
+
+  // Load needed material definition from media.geo file
+  create_materials_from_media_file();
+
+  // Get the GeoManager for later usage
+  gGeoMan = (TGeoManager*) gROOT->FindObject("FAIRGeom");
+  gGeoMan->SetVisLevel(5);  // 2 = super modules
+  gGeoMan->SetVisOption(0);
+
+  // Create the top volume
+  /*
+  TGeoBBox* topbox= new TGeoBBox("", 1000., 1000., 1000.);
+  TGeoVolume* top = new TGeoVolume("top", topbox, gGeoMan->GetMedium("air"));
+  gGeoMan->SetTopVolume(top);
+  */
+
+  TGeoVolume* top = new TGeoVolumeAssembly("TOP");
+  gGeoMan->SetTopVolume(top);
+
+  TGeoRotation* tof_rotation = new TGeoRotation();
+  tof_rotation->RotateY(0.);  // angle with respect to beam axis
+    //tof_rotation->RotateZ(   0 );   // electronics on  9 o'clock position = +x
+  //  tof_rotation->RotateZ(   0 );   // electronics on  9 o'clock position = +x
+  //  tof_rotation->RotateZ(  90 );   // electronics on 12 o'clock position (top)
+  //  tof_rotation->RotateZ( 180 );   // electronics on  3 o'clock position = -x
+  //  tof_rotation->RotateZ( 270 );   // electronics on  6 o'clock position (bottom)
+
+  TGeoVolume* tof = new TGeoVolumeAssembly(geoVersion);
+  //  top->AddNode(tof, 1, tof_rotation);
+  top->AddNode(tof, 1);
+
+  TGeoVolume* tofstand = new TGeoVolumeAssembly(geoVersionStand);
+  // Mar 2020 run
+  TGeoTranslation* stand_trans_local =
+    new TGeoTranslation("", TOF_X_Front_Stand, 0., 0.);
+  TGeoTranslation* stand_trans =
+    new TGeoTranslation("", 0., 0., TOF_Z_Front_Stand);    
+  // Nov 2019 run
+  // TGeoTranslation*  stand_trans   = new TGeoTranslation("", 12., 0., TOF_Z_Front_Stand);
+  // TGeoTranslation*  stand_trans   = new TGeoTranslation("",  0., 0., TOF_Z_Front_Stand);
+  TGeoRotation* stand_rot = new TGeoRotation();
+  stand_rot->RotateY(0.55);
+  TGeoCombiTrans* stand_combi_trans_local =
+    new TGeoCombiTrans(*stand_trans_local, *stand_rot);
+  TGeoCombiTrans* stand_combi_trans =
+    new TGeoCombiTrans(*stand_trans, *tof_rotation);    
+  //tof->AddNode(tofstand, 1, stand_combi_trans);
+  tof->AddNode(tofstand, 1, stand_combi_trans_local);
+  //tof->AddNode(tofstand, 1);
+
+  for (Int_t counterType = 0; counterType < NumberOfDifferentCounterTypes;
+       counterType++) {
+    gCounter[counterType] = create_new_counter(counterType);
+  }
+
+  for (Int_t moduleType = 0; moduleType < NofModuleTypes; moduleType++) {
+    gModules[moduleType] = create_new_tof_module(moduleType);
+    gModules[moduleType]->SetVisContainers(1);
+  }
+
+  // no pole
+  //  gPole = create_tof_pole();
+
+  //  position_side_tof_modules(1);  // keep order !!
+  //  position_inner_tof_modules(2);
+  position_inner_tof_modules(3);
+  position_Dia(1);
+  position_Star2(1);
+  //  position_cer_modules(3);
+  //  position_CERN(1);
+  position_Buc(1);
+
+  cout << "Outer Types " << Outer_Module_Types[0][0] << ", "
+       << Outer_Module_Types[1][0] << ", col=1:  " << Outer_Module_Types[0][1]
+       << ", " << Outer_Module_Types[1][1] << endl;
+  cout << "Outer Number " << Outer_Module_Number[0][0] << ", "
+       << Outer_Module_Number[1][0] << ", col=1:  " << Outer_Module_Number[0][1]
+       << ", " << Outer_Module_Number[1][1] << endl;
+  //  position_outer_tof_modules(4);
+  // position_tof_poles(0);
+  // position_tof_bars(0);
+
+  gGeoMan->CloseGeometry();
+  gGeoMan->CheckOverlaps(0.001);
+  gGeoMan->PrintOverlaps();
+  gGeoMan->CheckOverlaps(0.001, "s");
+  gGeoMan->PrintOverlaps();
+  gGeoMan->Test();
+
+  tof->Export(FileNameSim);
+  TFile* geoFile = new TFile(FileNameSim, "UPDATE");
+  stand_combi_trans->Write();
+  geoFile->Close();
+
+  /*
+  TFile* outfile1 = new TFile(FileNameSim,"RECREATE");
+  top->Write();
+  //gGeoMan->Write();
+  outfile1->Close();
+*/
+  //tof->RemoveNode((TGeoNode*)tofstand);
+  //top->AddNode(tof, 1, tof_rotation);
+  //tof->ReplaceNode((TGeoNode*)tofstand, 0, stand_combi_trans);
+  /*
+  CbmTransport run;
+  run.SetGeoFileName(FileNameGeo);
+  run.LoadSetup("setup_mcbm_tof_2020");
+  run.SetField(new CbmFieldConst());
+*/
+  //top->Export(FileNameGeo);
+
+  TFile* outfile2 = new TFile(FileNameGeo, "RECREATE");
+  gGeoMan->Write();
+  outfile2->Close();
+
+  dump_info_file();
+
+  top->SetVisContainers(1);
+  gGeoMan->SetVisLevel(5);
+  top->Draw("ogl");
+  //top->Draw();
+  //gModules[0]->Draw("ogl");
+  //  gModules[0]->Draw("");
+  gModules[0]->SetVisContainers(1);
+  //  gModules[1]->Draw("");
+  gModules[1]->SetVisContainers(1);
+  //gModules[5]->Draw("");
+  //  top->Raytrace();
+}
+
+void create_materials_from_media_file() {
+  // Use the FairRoot geometry interface to load the media which are already defined
+  FairGeoLoader* geoLoad    = new FairGeoLoader("TGeo", "FairGeoLoader");
+  FairGeoInterface* geoFace = geoLoad->getGeoInterface();
+  TString geoPath           = gSystem->Getenv("VMCWORKDIR");
+  TString geoFile           = geoPath + "/geometry/media.geo";
+  geoFace->setMediaFile(geoFile);
+  geoFace->readMedia();
+
+  // Read the required media and create them in the GeoManager
+  FairGeoMedia* geoMedia   = geoFace->getMedia();
+  FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();
+
+  FairGeoMedium* air          = geoMedia->getMedium("air");
+  FairGeoMedium* aluminium    = geoMedia->getMedium("aluminium");
+  FairGeoMedium* RPCgas       = geoMedia->getMedium("RPCgas");
+  FairGeoMedium* RPCgas_noact = geoMedia->getMedium("RPCgas_noact");
+  FairGeoMedium* RPCglass     = geoMedia->getMedium("RPCglass");
+  FairGeoMedium* carbon       = geoMedia->getMedium("carbon");
+
+  // include check if all media are found
+
+  geoBuild->createMedium(air);
+  geoBuild->createMedium(aluminium);
+  geoBuild->createMedium(RPCgas);
+  geoBuild->createMedium(RPCgas_noact);
+  geoBuild->createMedium(RPCglass);
+  geoBuild->createMedium(carbon);
+}
+
+TGeoVolume* create_counter(Int_t modType) {
+
+  //glass
+  Float_t gdx = Glass_X[modType];
+  Float_t gdy = Glass_Y[modType];
+  Float_t gdz = Glass_Z[modType];
+
+  //gas gap
+  Int_t nstrips = NumberOfReadoutStrips[modType];
+  Int_t ngaps   = NumberOfGaps[modType];
+
+
+  Float_t ggdx = GasGap_X[modType];
+  Float_t ggdy = GasGap_Y[modType];
+  Float_t ggdz = GasGap_Z[modType];
+  Float_t gsdx = ggdx / float(nstrips);
+
+  //single stack
+  Float_t dzpos     = gdz + ggdz;
+  Float_t startzpos = SingleStackStartPosition_Z[modType];
+
+  // electronics
+  //pcb dimensions
+  Float_t dxe  = Electronics_X[modType];
+  Float_t dye  = Electronics_Y[modType];
+  Float_t dze  = Electronics_Z[modType];
+  Float_t yele = (gdy + 0.1) / 2. + dye / 2.;
+
+  // needed materials
+  TGeoMedium* glassPlateVolMed  = gGeoMan->GetMedium(GlasMedium);
+  TGeoMedium* noActiveGasVolMed = gGeoMan->GetMedium(NoActivGasMedium);
+  TGeoMedium* activeGasVolMed   = gGeoMan->GetMedium(ActivGasMedium);
+  TGeoMedium* electronicsVolMed = gGeoMan->GetMedium(ElectronicsMedium);
+
+  // Single glass plate
+  TGeoBBox* glass_plate = new TGeoBBox("", gdx / 2., gdy / 2., gdz / 2.);
+  TGeoVolume* glass_plate_vol =
+    new TGeoVolume("tof_glass", glass_plate, glassPlateVolMed);
+  glass_plate_vol->SetLineColor(
+    kMagenta);                           // set line color for the glass plate
+  glass_plate_vol->SetTransparency(20);  // set transparency for the TOF
+  TGeoTranslation* glass_plate_trans = new TGeoTranslation("", 0., 0., 0.);
+
+  // Single gas gap
+  TGeoBBox* gas_gap = new TGeoBBox("", ggdx / 2., ggdy / 2., ggdz / 2.);
+  //TGeoVolume* gas_gap_vol =
+  //new TGeoVolume("tof_gas_gap", gas_gap, noActiveGasVolMed);
+  TGeoVolume* gas_gap_vol =
+    new TGeoVolume("tof_gas_active", gas_gap, activeGasVolMed);
+  gas_gap_vol->Divide("Strip", 1, nstrips, -ggdx / 2., 0);
+
+  gas_gap_vol->SetLineColor(kRed);   // set line color for the gas gap
+  gas_gap_vol->SetTransparency(70);  // set transparency for the TOF
+  TGeoTranslation* gas_gap_trans =
+    new TGeoTranslation("", 0., 0., (gdz + ggdz) / 2.);
+
+
+  // Single subdivided active gas gap
+  /*
+    TGeoBBox* gas_active = new TGeoBBox("", gsdx/2., ggdy/2., ggdz/2.);
+    TGeoVolume* gas_active_vol = 
+    new TGeoVolume("tof_gas_active", gas_active, activeGasVolMed);
+  gas_active_vol->SetLineColor(kBlack); // set line color for the gas gap
+  gas_active_vol->SetTransparency(70); // set transparency for the TOF
+  */
+
+  // Add glass plate, inactive gas gap and active gas gaps to a single stack
+  TGeoVolume* single_stack = new TGeoVolumeAssembly("single_stack");
+  single_stack->AddNode(glass_plate_vol, 0, glass_plate_trans);
+  single_stack->AddNode(gas_gap_vol, 0, gas_gap_trans);
+
+  /*
+  for (Int_t l=0; l<nstrips; l++){
+    TGeoTranslation* gas_active_trans 
+      = new TGeoTranslation("", -ggdx/2+(l+0.5)*gsdx, 0., 0.);
+    gas_gap_vol->AddNode(gas_active_vol, l, gas_active_trans);
+    //    single_stack->AddNode(gas_active_vol, l, gas_active_trans);
+  }
+  */
+
+  // Add 8 single stacks + one glass plate at the e09.750nd to a multi stack
+  TGeoVolume* multi_stack = new TGeoVolumeAssembly("multi_stack");
+  Int_t l;
+  for (l = 0; l < ngaps; l++) {
+    TGeoTranslation* single_stack_trans =
+      new TGeoTranslation("", 0., 0., startzpos + l * dzpos);
+    multi_stack->AddNode(single_stack, l, single_stack_trans);
+  }
+  TGeoTranslation* single_glass_back_trans =
+    new TGeoTranslation("", 0., 0., startzpos + ngaps * dzpos);
+  multi_stack->AddNode(glass_plate_vol, l, single_glass_back_trans);
+
+  // Add electronics above and below the glass stack to build a complete counter
+  TGeoVolume* counter                = new TGeoVolumeAssembly("counter");
+  TGeoTranslation* multi_stack_trans = new TGeoTranslation("", 0., 0., 0.);
+  counter->AddNode(multi_stack, l, multi_stack_trans);
+
+  TGeoBBox* pcb       = new TGeoBBox("", dxe / 2., dye / 2., dze / 2.);
+  TGeoVolume* pcb_vol = new TGeoVolume("pcb", pcb, electronicsVolMed);
+  pcb_vol->SetLineColor(kCyan);  // set line color for the gas gap
+  pcb_vol->SetTransparency(10);  // set transparency for the TOF
+  for (Int_t l = 0; l < 2; l++) {
+    yele *= -1.;
+    TGeoTranslation* pcb_trans = new TGeoTranslation("", 0., yele, 0.);
+    counter->AddNode(pcb_vol, l, pcb_trans);
+  }
+
+  return counter;
+}
+
+TGeoVolume* create_new_counter(Int_t modType) {
+
+  //glass
+  Float_t gdx = Glass_X[modType];
+  Float_t gdy = Glass_Y[modType];
+  Float_t gdz = Glass_Z[modType];
+
+  //gas gap
+  Int_t nstrips = NumberOfReadoutStrips[modType];
+  Int_t ngaps   = NumberOfGaps[modType];
+
+
+  Float_t ggdx = GasGap_X[modType];
+  Float_t ggdy = GasGap_Y[modType];
+  Float_t ggdz = GasGap_Z[modType];
+  Float_t gsdx = ggdx / (Float_t)(nstrips);
+
+  // electronics
+  //pcb dimensions
+  Float_t dxe  = Electronics_X[modType];
+  Float_t dye  = Electronics_Y[modType];
+  Float_t dze  = Electronics_Z[modType];
+  Float_t yele = gdy / 2. + dye / 2.;
+
+  // counter size (calculate from glas, gap and electronics sizes)
+  Float_t cdx = TMath::Max(gdx, ggdx);
+  cdx         = TMath::Max(cdx, dxe) + 0.2;
+  Float_t cdy = TMath::Max(gdy, ggdy) + 2 * dye + 0.2;
+  Float_t cdz = ngaps * ggdz + (ngaps + 1) * gdz
+                + 0.2;  // ngaps * (gdz+ggdz) + gdz + 0.2; // ok
+
+  //calculate thickness and first position in counter of single stack
+  Float_t dzpos = gdz + ggdz;
+  Float_t startzposglas =
+    -ngaps * (gdz + ggdz)
+    / 2.;  // -cdz/2.+0.1+gdz/2.; // ok  // (-cdz+gdz)/2.; // not ok
+  Float_t startzposgas =
+    startzposglas + gdz / 2. + ggdz / 2.;  // -cdz/2.+0.1+gdz   +ggdz/2.;  // ok
+
+
+  // needed materials
+  TGeoMedium* glassPlateVolMed  = gGeoMan->GetMedium(GlasMedium);
+  TGeoMedium* noActiveGasVolMed = gGeoMan->GetMedium(NoActivGasMedium);
+  TGeoMedium* activeGasVolMed   = gGeoMan->GetMedium(ActivGasMedium);
+  TGeoMedium* electronicsVolMed = gGeoMan->GetMedium(ElectronicsMedium);
+
+
+  // define counter volume
+  TGeoBBox* counter_box = new TGeoBBox("", cdx / 2., cdy / 2., cdz / 2.);
+  TGeoVolume* counter =
+    new TGeoVolume("counter", counter_box, noActiveGasVolMed);
+  counter->SetLineColor(kRed);   // set line color for the counter
+  counter->SetTransparency(70);  // set transparency for the TOF
+
+  // define single glass plate volume
+  TGeoBBox* glass_plate = new TGeoBBox("", gdx / 2., gdy / 2., gdz / 2.);
+  TGeoVolume* glass_plate_vol =
+    new TGeoVolume("tof_glass", glass_plate, glassPlateVolMed);
+  glass_plate_vol->SetLineColor(
+    kMagenta);                           // set line color for the glass plate
+  glass_plate_vol->SetTransparency(20);  // set transparency for the TOF
+  // define single gas gap volume
+  TGeoBBox* gas_gap       = new TGeoBBox("", ggdx / 2., ggdy / 2., ggdz / 2.);
+  TGeoVolume* gas_gap_vol = new TGeoVolume("Gap", gas_gap, activeGasVolMed);
+  gas_gap_vol->Divide("Cell", 1, nstrips, -ggdx / 2., 0);
+  gas_gap_vol->SetLineColor(kRed);   // set line color for the gas gap
+  gas_gap_vol->SetTransparency(99);  // set transparency for the TOF
+
+  // place 8 gas gaps and 9 glas plates in the counter
+  for (Int_t igap = 0; igap <= ngaps; igap++) {
+    // place (ngaps+1) glass plates
+    Float_t zpos_glas = startzposglas + igap * dzpos;
+    TGeoTranslation* glass_plate_trans =
+      new TGeoTranslation("", 0., 0., zpos_glas);
+    counter->AddNode(glass_plate_vol, igap, glass_plate_trans);
+    // place ngaps gas gaps
+    if (igap < ngaps) {
+      Float_t zpos_gas = startzposgas + igap * dzpos;
+      TGeoTranslation* gas_gap_trans =
+        new TGeoTranslation("", 0., 0., zpos_gas);
+      counter->AddNode(gas_gap_vol, igap, gas_gap_trans);
+    }
+    //    cout <<"Zpos(Glas): "<< zpos_glas << endl;
+    //    cout <<"Zpos(Gas): "<< zpos_gas << endl;
+  }
+
+  // create and place the electronics above and below the glas stack
+  TGeoBBox* pcb       = new TGeoBBox("", dxe / 2., dye / 2., dze / 2.);
+  TGeoVolume* pcb_vol = new TGeoVolume("pcb", pcb, electronicsVolMed);
+  pcb_vol->SetLineColor(kYellow);  // kCyan); // set line color for electronics
+  pcb_vol->SetTransparency(10);    // set transparency for the TOF
+  for (Int_t l = 0; l < 2; l++) {
+    yele *= -1.;
+    TGeoTranslation* pcb_trans = new TGeoTranslation("", 0., yele, 0.);
+    counter->AddNode(pcb_vol, l, pcb_trans);
+  }
+
+
+  return counter;
+}
+
+TGeoVolume* create_tof_module(Int_t modType) {
+  Int_t cType         = CounterTypeInModule[modType];
+  Float_t dx          = Module_Size_X[modType];
+  Float_t dy          = Module_Size_Y[modType];
+  Float_t dz          = Module_Size_Z[modType];
+  Float_t width_aluxl = Module_Thick_Alu_X_left;
+  Float_t width_aluxr = Module_Thick_Alu_X_right;
+  Float_t width_aluy  = Module_Thick_Alu_Y;
+  Float_t width_aluz  = Module_Thick_Alu_Z;
+
+  Float_t shift_gas_box =
+    (Module_Thick_Alu_X_right - Module_Thick_Alu_X_left) / 2;
+
+  Float_t dxpos     = CounterXDistance[modType];
+  Float_t startxpos = CounterXStartPosition[modType];
+  Float_t dzoff     = CounterZDistance[modType];
+  Float_t rotangle  = CounterRotationAngle[modType];
+
+  TGeoMedium* boxVolMed         = gGeoMan->GetMedium(BoxVolumeMedium);
+  TGeoMedium* noActiveGasVolMed = gGeoMan->GetMedium(NoActivGasMedium);
+
+  TString moduleName = Form("module_%d", modType);
+  TGeoVolume* module = new TGeoVolumeAssembly(moduleName);
+
+  TGeoBBox* alu_box       = new TGeoBBox("", dx / 2., dy / 2., dz / 2.);
+  TGeoVolume* alu_box_vol = new TGeoVolume("alu_box", alu_box, boxVolMed);
+  alu_box_vol->SetLineColor(kGreen);  // set line color for the alu box
+  alu_box_vol->SetTransparency(20);   // set transparency for the TOF
+  TGeoTranslation* alu_box_trans = new TGeoTranslation("", 0., 0., 0.);
+  module->AddNode(alu_box_vol, 0, alu_box_trans);
+
+  TGeoBBox* gas_box = new TGeoBBox("",
+                                   (dx - (width_aluxl + width_aluxr)) / 2.,
+                                   (dy - 2 * width_aluy) / 2.,
+                                   (dz - 2 * width_aluz) / 2.);
+  TGeoVolume* gas_box_vol =
+    new TGeoVolume("gas_box", gas_box, noActiveGasVolMed);
+  gas_box_vol->SetLineColor(kYellow);  // set line color for the gas box
+  gas_box_vol->SetTransparency(70);    // set transparency for the TOF
+  TGeoTranslation* gas_box_trans =
+    new TGeoTranslation("", shift_gas_box, 0., 0.);
+  alu_box_vol->AddNode(gas_box_vol, 0, gas_box_trans);
+
+  for (Int_t j = 0; j < 5; j++) {  //loop over counters (modules)
+    Float_t zpos;
+    if (0 == modType) {
+      zpos = dzoff *= -1;
+    } else {
+      zpos = 0.;
+    }
+    //cout << "counter z position " << zpos << endl;
+    TGeoTranslation* counter_trans =
+      new TGeoTranslation("", startxpos + j * dxpos, 0.0, zpos);
+
+    TGeoRotation* counter_rot = new TGeoRotation();
+    counter_rot->RotateY(rotangle);
+    TGeoCombiTrans* counter_combi_trans =
+      new TGeoCombiTrans(*counter_trans, *counter_rot);
+    gas_box_vol->AddNode(gCounter[cType], j, counter_combi_trans);
+  }
+
+  return module;
+}
+
+TGeoVolume* create_new_tof_module(Int_t modType) {
+  Int_t cType         = CounterTypeInModule[modType];
+  Float_t dx          = Module_Size_X[modType];
+  Float_t dy          = Module_Size_Y[modType];
+  Float_t dz          = Module_Size_Z[modType];
+  Float_t width_aluxl = Module_Thick_Alu_X_left;
+  Float_t width_aluxr = Module_Thick_Alu_X_right;
+  Float_t width_aluy  = Module_Thick_Alu_Y;
+  Float_t width_aluz  = Module_Thick_Alu_Z;
+
+  Float_t shift_gas_box =
+    (Module_Thick_Alu_X_right - Module_Thick_Alu_X_left) / 2;
+
+  Float_t dxpos     = CounterXDistance[modType];
+  Float_t startxpos = CounterXStartPosition[modType];
+  Float_t dypos     = CounterYDistance[modType];
+  Float_t startypos = CounterYStartPosition[modType];
+  Float_t dzoff     = CounterZDistance[modType];
+  Float_t rotangle  = CounterRotationAngle[modType];
+
+  TGeoMedium* boxVolMed         = gGeoMan->GetMedium(BoxVolumeMedium);
+  TGeoMedium* noActiveGasVolMed = gGeoMan->GetMedium(NoActivGasMedium);
+
+  TString moduleName = Form("module_%d", modType);
+
+  TGeoBBox* module_box = new TGeoBBox("", dx / 2., dy / 2., dz / 2.);
+  TGeoVolume* module   = new TGeoVolume(moduleName, module_box, boxVolMed);
+  module->SetLineColor(kGreen);  // set line color for the alu box
+  module->SetTransparency(20);   // set transparency for the TOF
+
+  TGeoBBox* gas_box = new TGeoBBox("",
+                                   (dx - (width_aluxl + width_aluxr)) / 2.,
+                                   (dy - 2 * width_aluy) / 2.,
+                                   (dz - 2 * width_aluz) / 2.);
+  TGeoVolume* gas_box_vol =
+    new TGeoVolume("gas_box", gas_box, noActiveGasVolMed);
+  gas_box_vol->SetLineColor(kBlue);  // set line color for the alu box
+  gas_box_vol->SetTransparency(50);  // set transparency for the TOF
+  TGeoTranslation* gas_box_trans =
+    new TGeoTranslation("", shift_gas_box, 0., 0.);
+  module->AddNode(gas_box_vol, 0, gas_box_trans);
+
+  for (Int_t j = 0; j < NCounterInModule[modType];
+       j++) {  //loop over counters (modules)
+               //for (Int_t j=0; j< 1; j++){ //loop over counters (modules)
+    Float_t xpos, ypos, zpos;
+    if (0 == modType || 3 == modType || 4 == modType || 5 == modType) {
+      zpos = dzoff *= -1;
+    } else {
+      zpos = CounterZStartPosition[modType] + j * dzoff;
+    }
+    //cout << "counter z position " << zpos << endl;
+    xpos = startxpos + j * dxpos;
+    ypos = startypos + j * dypos;
+
+    TGeoTranslation* counter_trans = new TGeoTranslation("", xpos, ypos, zpos);
+
+    TGeoRotation* counter_rot = new TGeoRotation();
+    counter_rot->RotateY(rotangle);
+    TGeoCombiTrans* counter_combi_trans =
+      new TGeoCombiTrans(*counter_trans, *counter_rot);
+    gas_box_vol->AddNode(gCounter[cType], j, counter_combi_trans);
+  }
+
+  return module;
+}
+
+
+TGeoVolume* create_tof_pole() {
+  // needed materials
+  TGeoMedium* boxVolMed = gGeoMan->GetMedium(BoxVolumeMedium);
+  TGeoMedium* airVolMed = gGeoMan->GetMedium(KeepingVolumeMedium);
+
+  Float_t dx         = Pole_Size_X;
+  Float_t dy         = Pole_Size_Y;
+  Float_t dz         = Pole_Size_Z;
+  Float_t width_alux = Pole_Thick_X;
+  Float_t width_aluy = Pole_Thick_Y;
+  Float_t width_aluz = Pole_Thick_Z;
+
+  TGeoVolume* pole       = new TGeoVolumeAssembly("Pole");
+  TGeoBBox* pole_alu_box = new TGeoBBox("", dx / 2., dy / 2., dz / 2.);
+  TGeoVolume* pole_alu_vol =
+    new TGeoVolume("pole_alu", pole_alu_box, boxVolMed);
+  pole_alu_vol->SetLineColor(kGreen);  // set line color for the alu box
+  pole_alu_vol->SetTransparency(20);   // set transparency for the TOF
+  TGeoTranslation* pole_alu_trans = new TGeoTranslation("", 0., 0., 0.);
+  pole->AddNode(pole_alu_vol, 0, pole_alu_trans);
+
+  Float_t air_dx = dx / 2. - width_alux;
+  Float_t air_dy = dy / 2. - width_aluy;
+  Float_t air_dz = dz / 2. - width_aluz;
+
+  //  cout << "My pole." << endl;
+  if (air_dx <= 0.)
+    cout << "ERROR - No air volume in pole X, size: " << air_dx << endl;
+  if (air_dy <= 0.)
+    cout << "ERROR - No air volume in pole Y, size: " << air_dy << endl;
+  if (air_dz <= 0.)
+    cout << "ERROR - No air volume in pole Z, size: " << air_dz << endl;
+
+  if ((air_dx > 0.) && (air_dy > 0.)
+      && (air_dz > 0.))  // crate air volume only, if larger than zero
+  {
+    TGeoBBox* pole_air_box = new TGeoBBox("", air_dx, air_dy, air_dz);
+    //  TGeoBBox* pole_air_box = new TGeoBBox("", dx/2.-width_alux, dy/2.-width_aluy, dz/2.-width_aluz);
+    TGeoVolume* pole_air_vol =
+      new TGeoVolume("pole_air", pole_air_box, airVolMed);
+    pole_air_vol->SetLineColor(kYellow);  // set line color for the alu box
+    pole_air_vol->SetTransparency(70);    // set transparency for the TOF
+    TGeoTranslation* pole_air_trans = new TGeoTranslation("", 0., 0., 0.);
+    pole_alu_vol->AddNode(pole_air_vol, 0, pole_air_trans);
+  } else
+    cout << "Skipping pole_air_vol, no thickness: " << air_dx << " " << air_dy
+         << " " << air_dz << endl;
+
+  return pole;
+}
+
+TGeoVolume* create_tof_bar(Float_t dx, Float_t dy, Float_t dz) {
+  // needed materials
+  TGeoMedium* boxVolMed = gGeoMan->GetMedium(BoxVolumeMedium);
+  TGeoMedium* airVolMed = gGeoMan->GetMedium(KeepingVolumeMedium);
+
+  Float_t width_alux = Pole_Thick_X;
+  Float_t width_aluy = Pole_Thick_Y;
+  Float_t width_aluz = Pole_Thick_Z;
+
+  TGeoVolume* bar         = new TGeoVolumeAssembly("Bar");
+  TGeoBBox* bar_alu_box   = new TGeoBBox("", dx / 2., dy / 2., dz / 2.);
+  TGeoVolume* bar_alu_vol = new TGeoVolume("bar_alu", bar_alu_box, boxVolMed);
+  bar_alu_vol->SetLineColor(kGreen);  // set line color for the alu box
+  bar_alu_vol->SetTransparency(20);   // set transparency for the TOF
+  TGeoTranslation* bar_alu_trans = new TGeoTranslation("", 0., 0., 0.);
+  bar->AddNode(bar_alu_vol, 0, bar_alu_trans);
+
+  TGeoBBox* bar_air_box = new TGeoBBox(
+    "", dx / 2. - width_alux, dy / 2. - width_aluy, dz / 2. - width_aluz);
+  TGeoVolume* bar_air_vol = new TGeoVolume("bar_air", bar_air_box, airVolMed);
+  bar_air_vol->SetLineColor(kYellow);  // set line color for the alu box
+  bar_air_vol->SetTransparency(70);    // set transparency for the TOF
+  TGeoTranslation* bar_air_trans = new TGeoTranslation("", 0., 0., 0.);
+  bar_alu_vol->AddNode(bar_air_vol, 0, bar_air_trans);
+
+  return bar;
+}
+
+void position_tof_poles(Int_t modType) {
+
+  TGeoTranslation* pole_trans = NULL;
+
+  Int_t numPoles = 0;
+  for (Int_t i = 0; i < NumberOfPoles; i++) {
+    if (i < 2) {
+      pole_trans =
+        new TGeoTranslation("", -Pole_Offset + 2.0, 0., Pole_ZPos[i]);
+      gGeoMan->GetVolume(geoVersionStand)->AddNode(gPole, numPoles, pole_trans);
+      numPoles++;
+    } else {
+      Float_t xPos = Pole_Offset + Pole_Size_X / 2. + Pole_Col[i] * DxColl;
+      Float_t zPos = Pole_ZPos[i];
+      pole_trans   = new TGeoTranslation("", xPos, 0., zPos);
+      gGeoMan->GetVolume(geoVersionStand)->AddNode(gPole, numPoles, pole_trans);
+      numPoles++;
+
+      pole_trans = new TGeoTranslation("", -xPos, 0., zPos);
+      gGeoMan->GetVolume(geoVersionStand)->AddNode(gPole, numPoles, pole_trans);
+      numPoles++;
+    }
+    cout << " Position Pole " << numPoles << " at z=" << Pole_ZPos[i] << endl;
+  }
+}
+
+void position_tof_bars(Int_t modType) {
+
+  TGeoTranslation* bar_trans = NULL;
+
+  Int_t numBars = 0;
+  Int_t i;
+  Float_t xPos;
+  Float_t yPos;
+  Float_t zPos;
+
+  for (i = 0; i < NumberOfBars; i++) {
+
+    xPos = Bar_XPos[i];
+    zPos = Bar_ZPos[i];
+    yPos = Pole_Size_Y / 2. + Bar_Size_Y / 2.;
+
+    bar_trans = new TGeoTranslation("", xPos, yPos, zPos);
+    gGeoMan->GetVolume(geoVersionStand)->AddNode(gBar[i], numBars, bar_trans);
+    numBars++;
+
+    bar_trans = new TGeoTranslation("", xPos, -yPos, zPos);
+    gGeoMan->GetVolume(geoVersionStand)->AddNode(gBar[i], numBars, bar_trans);
+    numBars++;
+
+    bar_trans = new TGeoTranslation("", -xPos, yPos, zPos);
+    gGeoMan->GetVolume(geoVersionStand)->AddNode(gBar[i], numBars, bar_trans);
+    numBars++;
+
+    bar_trans = new TGeoTranslation("", -xPos, -yPos, zPos);
+    gGeoMan->GetVolume(geoVersionStand)->AddNode(gBar[i], numBars, bar_trans);
+    numBars++;
+  }
+  cout << " Position Bar " << numBars << " at z=" << Bar_ZPos[i] << endl;
+
+  // horizontal frame bars
+  i = NumberOfBars;
+  NumberOfBars++;
+  // no bar
+  //   gBar[i]=create_tof_bar(2.*xPos+Pole_Size_X,Bar_Size_Y,Bar_Size_Y);
+
+  zPos      = Pole_ZPos[0] + Pole_Size_Z / 2.;
+  bar_trans = new TGeoTranslation("", 0., yPos, zPos);
+  gGeoMan->GetVolume(geoVersionStand)->AddNode(gBar[i], numBars, bar_trans);
+  numBars++;
+
+  bar_trans = new TGeoTranslation("", 0., -yPos, zPos);
+  gGeoMan->GetVolume(geoVersionStand)->AddNode(gBar[i], numBars, bar_trans);
+  numBars++;
+}
+
+void position_inner_tof_modules(Int_t modNType) {
+  TGeoTranslation* module_trans = NULL;
+
+  //  Int_t numModules=(Int_t)( (Inner_Module_Last_Y_Position-Inner_Module_First_Y_Position)/Module_Size_Y[modType])+1;
+  Float_t yPos = Inner_Module_First_Y_Position;
+  Int_t ii     = 0;
+  Float_t xPos = Inner_Module_X_Offset;
+  Float_t zPos = Wall_Z_Position;
+
+  Pole_ZPos[NumberOfPoles] = zPos;
+  Pole_Col[NumberOfPoles]  = 0;
+  NumberOfPoles++;
+
+  Float_t DzPos = 0.;
+  for (Int_t j = 0; j < modNType; j++) {
+    if (Module_Size_Z[j] > DzPos) { DzPos = Module_Size_Z[j]; }
+  }
+  Pole_ZPos[NumberOfPoles] = zPos + DzPos;
+  Pole_Col[NumberOfPoles]  = 0;
+  NumberOfPoles++;
+
+  // for (Int_t j=0; j<modNType; j++){
+  // for (Int_t j=1; j<modNType; j++){
+  Int_t modType;
+  Int_t modNum;
+  for (Int_t j = 2; j < modNType;
+       j++) {  // place only M4 type modules (modNType == 2)
+               //DEDE
+    modType = Inner_Module_Types[j];
+    modNum  = 0;
+    //  for(Int_t i=0; i<Inner_Module_Number[j]; i++) {
+    //  for(Int_t i=0; i<1; i++) { // place 1x2 modules in the top and same in the bottom
+    for (Int_t i = 0; i < 2;
+         i++) {  // place 2x2 modules in the top and same in the bottom
+      ii++;
+      cout << "Inner ii " << ii << " Last " << Last_Size_Y << ", "
+           << Last_Over_Y << endl;
+      Float_t DeltaY = Module_Size_Y[modType] + Last_Size_Y
+                       - 2. * (Module_Over_Y[modType] + Last_Over_Y);
+      //    DeltaY = 1.5;
+      cout << "DeltaY " << DeltaY << endl;
+      yPos += DeltaY;
+      Last_Size_Y = Module_Size_Y[modType];
+      Last_Over_Y = Module_Over_Y[modType];
+      cout << "Position Inner Module " << i << " of " << Inner_Module_Number[j]
+           << " Type " << modType << " at Y = " << yPos
+           << " Ysize = " << Module_Size_Y[modType] << " DeltaY = " << DeltaY
+           << endl;
+
+      ///    module_trans = new TGeoTranslation("", xPos, yPos, zPos);
+      ///    gGeoMan->GetVolume(geoVersionStand)->AddNode(gModules[modType], modNum, module_trans);
+      ///    modNum++;
+      ///    module_trans = new TGeoTranslation("", xPos, -yPos, zPos);
+      ///    gGeoMan->GetVolume(geoVersionStand)->AddNode(gModules[modType], modNum, module_trans);
+      ///    modNum++;
+      //    //    if (ii>0) {
+      //    if (ii>1) {
+      //      module_trans
+      //	= new TGeoTranslation("", xPos, yPos-DeltaY/2, zPos+Module_Size_Z[modType]);
+      //      gGeoMan->GetVolume(geoVersionStand)->AddNode(gModules[modType], modNum, module_trans);
+      //    modNum++;
+      //      module_trans
+      //	= new TGeoTranslation("", xPos, -(yPos-DeltaY/2), zPos+Module_Size_Z[modType]);
+      //      gGeoMan->GetVolume(geoVersionStand)->AddNode(gModules[modType], modNum, module_trans);
+      //    modNum++;
+      //    }
+    }
+  }
+  // module_trans = new TGeoTranslation("", xPos, -49-3, zPos);
+
+  // Mar2019 setup
+  const Int_t NModules           = 5;
+  xPos                           = 0.;
+  yPos                           = 0.;
+  zPos                           = TOF_Z_Front;
+  const Double_t ModDx[NModules] = {0., 0., 0., 49.8, 49.8};
+  //const Double_t ModDx[NModules]= { 1.5,    0., -1.5, 49.8, 55.8};
+  const Double_t ModDy[NModules]     = {0., 0., 0., 0., 0.};
+  const Double_t ModDz[NModules]     = {0., 16.5, 34., 0., 34.};
+  const Double_t ModAng[NModules]    = {-90., -90., -90., -90., -90.0};
+  TGeoRotation* module_rot           = NULL;
+  TGeoCombiTrans* module_combi_trans = NULL;
+
+  for (Int_t iMod = 0; iMod < NModules; iMod++) {
+    module_trans = new TGeoTranslation(
+      "", xPos + ModDx[iMod], yPos + ModDy[iMod], zPos + ModDz[iMod]);
+    module_rot = new TGeoRotation();
+    module_rot->RotateZ(ModAng[iMod]);
+    module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+    gGeoMan->GetVolume(geoVersionStand)
+      ->AddNode(gModules[modType], modNum, module_combi_trans);
+    modNum++;
+  }
+
+
+  /*
+ module_trans = new TGeoTranslation("", xPos, 0, zPos+16.5);
+ gGeoMan->GetVolume(geoVersionStand)->AddNode(gModules[modType], modNum, module_trans);
+ modNum++;
+ 
+ // module_trans = new TGeoTranslation("", xPos, 49+3, zPos);
+ module_trans = new TGeoTranslation("", xPos, 0, zPos+16.5+17.5);
+ gGeoMan->GetVolume(geoVersionStand)->AddNode(gModules[modType], modNum, module_trans);
+ modNum++;
+
+ // module_trans = new TGeoTranslation("", xPos,-26, zPos+Module_Size_Z[modType]);
+ module_trans = new TGeoTranslation("", xPos, -49.8, zPos);
+ gGeoMan->GetVolume(geoVersionStand)->AddNode(gModules[modType], modNum, module_trans);
+ modNum++;
+
+ // module_trans = new TGeoTranslation("", xPos, 26, zPos+Module_Size_Z[modType]);
+ module_trans = new TGeoTranslation("", xPos, -49.8, zPos+16.5);
+ gGeoMan->GetVolume(geoVersionStand)->AddNode(gModules[modType], modNum, module_trans);
+ modNum++;
+ */
+}
+
+
+void position_Dia(Int_t modNType) {
+  TGeoTranslation* module_trans = NULL;
+  TGeoRotation* module_rot      = new TGeoRotation();
+  module_rot->RotateZ(Dia_rotate_Z);
+  TGeoCombiTrans* module_combi_trans = NULL;
+
+  //  Int_t numModules=(Int_t)( (Inner_Module_Last_Y_Position-Inner_Module_First_Y_Position)/Module_Size_Y[modType])+1;
+  Float_t yPos = Dia_First_Y_Position;
+  Int_t ii     = 0;
+  Float_t xPos = Dia_X_Offset;
+  Float_t zPos = Dia_Z_Position;
+
+  Int_t modNum = 0;
+  for (Int_t j = 0; j < modNType; j++) {
+    Int_t modType = Dia_Types[j];
+    for (Int_t i = 0; i < Dia_Number[j]; i++) {
+      ii++;
+      module_trans       = new TGeoTranslation("", xPos, yPos, zPos);
+      module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+      gGeoMan->GetVolume(geoVersionStand)
+        ->AddNode(gModules[modType], modNum, module_combi_trans);
+      modNum++;
+    }
+  }
+}
+
+void position_Star2(Int_t modNType) {
+  TGeoTranslation* module_trans = NULL;
+  TGeoRotation* module_rot      = new TGeoRotation();
+  module_rot->RotateZ(Star2_rotate_Z);
+  TGeoCombiTrans* module_combi_trans = NULL;
+
+  Float_t yPos = Star2_First_Y_Position;
+  Float_t zPos = Star2_First_Z_Position;
+  Int_t ii     = 0;
+
+  Int_t modNum = 0;
+  for (Int_t j = 0; j < modNType; j++) {
+    Int_t modType = Star2_Types[j];
+    Float_t xPos  = Star2_X_Offset[j];
+    for (Int_t i = 0; i < Star2_Number[j]; i++) {
+      ii++;
+      module_trans       = new TGeoTranslation("", xPos, yPos, zPos);
+      module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+      gGeoMan->GetVolume(geoVersionStand)
+        ->AddNode(gModules[modType], modNum, module_combi_trans);
+      modNum++;
+      yPos += Star2_Delta_Y_Position;
+      zPos += Star2_Delta_Z_Position;
+    }
+  }
+}
+
+void position_Buc(Int_t modNType) {
+  TGeoTranslation* module_trans = NULL;
+  TGeoRotation* module_rot      = new TGeoRotation();
+  module_rot->RotateZ(Buc_rotate_Z);
+  TGeoCombiTrans* module_combi_trans = NULL;
+
+  Float_t yPos = Buc_First_Y_Position;
+  Float_t zPos = Buc_First_Z_Position;
+  Int_t ii     = 0;
+
+  Int_t modNum = 0;
+  for (Int_t j = 0; j < modNType; j++) {
+    Int_t modType = Buc_Types[j];
+    Float_t xPos  = Buc_X_Offset[j];
+    for (Int_t i = 0; i < Buc_Number[j]; i++) {
+      ii++;
+      module_trans       = new TGeoTranslation("", xPos, yPos, zPos);
+      module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+      gGeoMan->GetVolume(geoVersionStand)
+        ->AddNode(gModules[modType], modNum, module_combi_trans);
+      modNum++;
+      yPos += Buc_Delta_Y_Position;
+      zPos += Buc_Delta_Z_Position;
+    }
+  }
+}
+
+void position_cer_modules(Int_t modNType) {
+  Int_t ii     = 0;
+  Int_t modNum = 0;
+  for (Int_t j = 1; j < modNType; j++) {
+    Int_t modType                 = Cer_Types[j];
+    Float_t xPos                  = Cer_X_Position[j];
+    Float_t yPos                  = Cer_Y_Position[j];
+    Float_t zPos                  = Cer_Z_Position[j];
+    TGeoTranslation* module_trans = NULL;
+    TGeoRotation* module_rot =
+      new TGeoRotation(Form("Cer%d", j), Cer_rotate_Z[j], -MeanTheta, 0.);
+    // module_rot->RotateZ(Cer_rotate_Z[j]);
+    TGeoCombiTrans* module_combi_trans = NULL;
+
+    for (Int_t i = 0; i < Cer_Number[j]; i++) {
+      ii++;
+      cout << "Position Ceramic Module " << i << " of " << Cer_Number[j]
+           << " Type " << modType << " at X = " << xPos << ", Y = " << yPos
+           << ", Z = " << zPos << endl;
+      // Front staggered module (Top if pair), top
+      module_trans       = new TGeoTranslation("", xPos, yPos, zPos);
+      module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+      gGeoMan->GetVolume(geoVersionStand)
+        ->AddNode(gModules[modType], modNum, module_combi_trans);
+      //    modNum++;
+    }
+  }
+}
+
+void position_CERN(Int_t modNType) {
+  TGeoTranslation* module_trans = NULL;
+  TGeoRotation* module_rot      = new TGeoRotation();
+  module_rot->RotateZ(CERN_rotate_Z);
+  TGeoCombiTrans* module_combi_trans = NULL;
+
+  //  Int_t numModules=(Int_t)( (Inner_Module_Last_Y_Position-Inner_Module_First_Y_Position)/Module_Size_Y[modType])+1;
+  Float_t yPos = CERN_First_Y_Position;
+  Int_t ii     = 0;
+  Float_t xPos = CERN_X_Offset;
+  Float_t zPos = CERN_Z_Position;
+
+  for (Int_t j = 0; j < modNType; j++) {
+    Int_t modType = CERN_Types[j];
+    Int_t modNum  = 0;
+    for (Int_t i = 0; i < CERN_Number[j]; i++) {
+      ii++;
+      module_trans       = new TGeoTranslation("", xPos, yPos, zPos);
+      module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+      gGeoMan->GetVolume(geoVersionStand)
+        ->AddNode(gModules[modType], modNum, module_combi_trans);
+      modNum++;
+    }
+  }
+}
+
+void position_side_tof_modules(Int_t modNType) {
+  TGeoTranslation* module_trans = NULL;
+  TGeoRotation* module_rot      = new TGeoRotation();
+  module_rot->RotateZ(180.);
+  TGeoCombiTrans* module_combi_trans = NULL;
+
+  //  Int_t numModules=(Int_t)( (Inner_Module_Last_Y_Position-Inner_Module_First_Y_Position)/Module_Size_Y[modType])+1;
+  Float_t yPos = 0.;  //Inner_Module_First_Y_Position;
+  Int_t ii     = 0;
+  for (Int_t j = 0; j < modNType; j++) {
+    Int_t modType = InnerSide_Module_Types[j];
+    Int_t modNum  = 0;
+    for (Int_t i = 0; i < InnerSide_Module_Number[j]; i++) {
+      ii++;
+      cout << "InnerSide ii " << ii << " Last " << Last_Size_Y << ","
+           << Last_Over_Y << endl;
+      Float_t DeltaY = Module_Size_Y[modType] + Last_Size_Y
+                       - 2. * (Module_Over_Y[modType] + Last_Over_Y);
+      if (ii > 1) { yPos += DeltaY; }
+      Last_Size_Y  = Module_Size_Y[modType];
+      Last_Over_Y  = Module_Over_Y[modType];
+      Float_t xPos = InnerSide_Module_X_Offset;
+      Float_t zPos = Wall_Z_Position;
+      cout << "Position InnerSide Module " << i << " of "
+           << InnerSide_Module_Number[j] << " Type " << modType
+           << " at Y = " << yPos << " Ysize = " << Module_Size_Y[modType]
+           << " DeltaY = " << DeltaY << endl;
+
+      module_trans = new TGeoTranslation("", xPos, yPos, zPos);
+      gGeoMan->GetVolume(geoVersionStand)
+        ->AddNode(gModules[modType], modNum, module_trans);
+      modNum++;
+
+      module_trans       = new TGeoTranslation("", -xPos, yPos, zPos);
+      module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+      gGeoMan->GetVolume(geoVersionStand)
+        ->AddNode(gModules[modType], modNum, module_combi_trans);
+      modNum++;
+
+      if (ii > 1) {
+        module_trans = new TGeoTranslation("", xPos, -yPos, zPos);
+        gGeoMan->GetVolume(geoVersionStand)
+          ->AddNode(gModules[modType], modNum, module_trans);
+        modNum++;
+
+        module_trans       = new TGeoTranslation("", -xPos, -yPos, zPos);
+        module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+        gGeoMan->GetVolume(geoVersionStand)
+          ->AddNode(gModules[modType], modNum, module_combi_trans);
+        modNum++;
+
+        module_trans = new TGeoTranslation(
+          "", xPos, yPos - DeltaY / 2, zPos + Module_Size_Z[modType]);
+        gGeoMan->GetVolume(geoVersionStand)
+          ->AddNode(gModules[modType], modNum, module_trans);
+        modNum++;
+
+        module_trans = new TGeoTranslation(
+          "", -xPos, yPos - DeltaY / 2, zPos + Module_Size_Z[modType]);
+        module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+        gGeoMan->GetVolume(geoVersionStand)
+          ->AddNode(gModules[modType], modNum, module_combi_trans);
+        modNum++;
+
+        module_trans = new TGeoTranslation(
+          "", xPos, -(yPos - DeltaY / 2), zPos + Module_Size_Z[modType]);
+        gGeoMan->GetVolume(geoVersionStand)
+          ->AddNode(gModules[modType], modNum, module_trans);
+        modNum++;
+
+        module_trans = new TGeoTranslation(
+          "", -xPos, -(yPos - DeltaY / 2), zPos + Module_Size_Z[modType]);
+        module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+        gGeoMan->GetVolume(geoVersionStand)
+          ->AddNode(gModules[modType], modNum, module_combi_trans);
+        modNum++;
+      }
+    }
+  }
+}
+
+void position_outer_tof_modules(Int_t nCol)  //modType, Int_t col1, Int_t col2)
+{
+  TGeoTranslation* module_trans = NULL;
+  TGeoRotation* module_rot      = new TGeoRotation();
+  module_rot->RotateZ(180.);
+  TGeoCombiTrans* module_combi_trans = NULL;
+
+  //  Int_t numModules=(Int_t)( (Outer_Module_Last_Y_Position-Outer_Module_First_Y_Position)/Module_Size_Y[modType])+1;
+
+  Int_t modNum[NofModuleTypes];
+  for (Int_t k = 0; k < NofModuleTypes; k++) {
+    modNum[k] = 0;
+  }
+
+  Float_t zPos = Wall_Z_Position;
+  for (Int_t j = 0; j < nCol; j++) {
+    Float_t xPos  = Outer_Module_X_Offset + ((j + 1) * DxColl);
+    Last_Size_Y   = 0.;
+    Last_Over_Y   = 0.;
+    Float_t yPos  = 0.;
+    Int_t ii      = 0;
+    Float_t DzPos = 0.;
+    for (Int_t k = 0; k < Outer_Module_NTypes; k++) {
+      Int_t modType = Outer_Module_Types[k][j];
+      if (Module_Size_Z[modType] > DzPos) {
+        if (Outer_Module_Number[k][j] > 0) { DzPos = Module_Size_Z[modType]; }
+      }
+    }
+
+    zPos -= 2. * DzPos;  //((j+1)*2*Module_Size_Z[modType]);
+
+    Pole_ZPos[NumberOfPoles] = zPos;
+    Pole_Col[NumberOfPoles]  = j + 1;
+    NumberOfPoles++;
+    Pole_ZPos[NumberOfPoles] = zPos + DzPos;
+    Pole_Col[NumberOfPoles]  = j + 1;
+    NumberOfPoles++;
+    //if (j+1==nCol) {
+    if (1) {
+      Pole_ZPos[NumberOfPoles] = Pole_ZPos[0];
+      Pole_Col[NumberOfPoles]  = j + 1;
+      NumberOfPoles++;
+
+      Bar_Size_Z         = Pole_ZPos[0] - zPos;
+      gBar[NumberOfBars] = create_tof_bar(Bar_Size_X, Bar_Size_Y, Bar_Size_Z);
+      Bar_ZPos[NumberOfBars] = zPos + Bar_Size_Z / 2. - Pole_Size_Z / 2.;
+      Bar_XPos[NumberOfBars] = xPos + Pole_Offset;
+      NumberOfBars++;
+    }
+
+    for (Int_t k = 0; k < Outer_Module_NTypes; k++) {
+      Int_t modType    = Outer_Module_Types[k][j];
+      Int_t numModules = Outer_Module_Number[k][j];
+
+      cout << " Outer: position " << numModules << " of type " << modType
+           << " in col " << j << " at z = " << zPos << ", DzPos = " << DzPos
+           << endl;
+      for (Int_t i = 0; i < numModules; i++) {
+        ii++;
+        cout << "Outer ii " << ii << " Last " << Last_Size_Y << ","
+             << Last_Over_Y << endl;
+        Float_t DeltaY = Module_Size_Y[modType] + Last_Size_Y
+                         - 2. * (Module_Over_Y[modType] + Last_Over_Y);
+        if (ii > 1) { yPos += DeltaY; }
+        Last_Size_Y = Module_Size_Y[modType];
+        Last_Over_Y = Module_Over_Y[modType];
+        cout << "Position Outer Module " << i << " of "
+             << Outer_Module_Number[k][j] << " Type " << modType << "(#"
+             << modNum[modType] << ") "
+             << " at Y = " << yPos << " Ysize = " << Module_Size_Y[modType]
+             << " DeltaY = " << DeltaY << endl;
+
+        module_trans = new TGeoTranslation("", xPos, yPos, zPos);
+        gGeoMan->GetVolume(geoVersionStand)
+          ->AddNode(gModules[modType], modNum[modType], module_trans);
+        modNum[modType]++;
+
+        module_trans       = new TGeoTranslation("", -xPos, yPos, zPos);
+        module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+        gGeoMan->GetVolume(geoVersionStand)
+          ->AddNode(gModules[modType], modNum[modType], module_combi_trans);
+        modNum[modType]++;
+
+        if (ii > 1) {
+          module_trans = new TGeoTranslation("", xPos, -yPos, zPos);
+          gGeoMan->GetVolume(geoVersionStand)
+            ->AddNode(gModules[modType], modNum[modType], module_trans);
+          modNum[modType]++;
+          module_trans       = new TGeoTranslation("", -xPos, -yPos, zPos);
+          module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+          gGeoMan->GetVolume(geoVersionStand)
+            ->AddNode(gModules[modType], modNum[modType], module_combi_trans);
+          modNum[modType]++;
+
+          // second layer
+          module_trans =
+            new TGeoTranslation("", xPos, yPos - DeltaY / 2., zPos + DzPos);
+          gGeoMan->GetVolume(geoVersionStand)
+            ->AddNode(gModules[modType], modNum[modType], module_trans);
+          modNum[modType]++;
+          module_trans =
+            new TGeoTranslation("", -xPos, yPos - DeltaY / 2., zPos + DzPos);
+          module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+          gGeoMan->GetVolume(geoVersionStand)
+            ->AddNode(gModules[modType], modNum[modType], module_combi_trans);
+          modNum[modType]++;
+
+          module_trans =
+            new TGeoTranslation("", xPos, -(yPos - DeltaY / 2.), zPos + DzPos);
+          gGeoMan->GetVolume(geoVersionStand)
+            ->AddNode(gModules[modType], modNum[modType], module_trans);
+          modNum[modType]++;
+          module_trans =
+            new TGeoTranslation("", -xPos, -(yPos - DeltaY / 2.), zPos + DzPos);
+          module_combi_trans = new TGeoCombiTrans(*module_trans, *module_rot);
+          gGeoMan->GetVolume(geoVersionStand)
+            ->AddNode(gModules[modType], modNum[modType], module_combi_trans);
+          modNum[modType]++;
+        }
+      }
+    }
+  }
+}
+
+
+void dump_info_file() {
+  TDatime datetime;  // used to get timestamp
+
+  printf("writing info file: %s\n", FileNameInfo.Data());
+
+  FILE* ifile;
+  ifile = fopen(FileNameInfo.Data(), "w");
+
+  if (ifile == NULL) {
+    printf("error opening %s\n", FileNameInfo.Data());
+    exit(1);
+  }
+
+  fprintf(ifile, "#\n##   %s information file\n#\n\n", geoVersion.Data());
+
+  fprintf(ifile, "# created %d\n\n", datetime.GetDate());
+
+  fprintf(ifile, "# TOF setup\n");
+  if (TOF_Z_Front == 450) fprintf(ifile, "SIS 100 hadron setup\n");
+  if (TOF_Z_Front == 600) fprintf(ifile, "SIS 100 electron\n");
+  if (TOF_Z_Front == 650) fprintf(ifile, "SIS 100 muon\n");
+  if (TOF_Z_Front == 880) fprintf(ifile, "SIS 300 electron\n");
+  if (TOF_Z_Front == 1020) fprintf(ifile, "SIS 300 muon\n");
+  fprintf(ifile, "\n");
+
+  const Float_t TOF_Z_Back =
+    Wall_Z_Position + 1.5 * Module_Size_Z[0];  // back of TOF wall
+
+  fprintf(ifile, "# envelope\n");
+  // Show extension of TRD
+  fprintf(ifile, "%7.2f cm   start of TOF (z)\n", TOF_Z_Front);
+  fprintf(ifile, "%7.2f cm   end   of TOF (z)\n", TOF_Z_Back);
+  fprintf(ifile, "\n");
+
+  // Layer thickness
+  fprintf(ifile, "# central tower position\n");
+  fprintf(ifile,
+          "%7.2f cm   center of staggered, front RPC cell at x=0\n",
+          Wall_Z_Position);
+  fprintf(ifile, "\n");
+
+  fclose(ifile);
+}