diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 82c322f91f04b096d33d0e2d067837297594146a..164594d57926e8e4124046455457f139a4d8db64 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -47,9 +47,20 @@ common:detector:
     - mkdir -p config
     - print_env.sh
 
+.phy_benchmark:
+  needs:
+    - ["common:detector"]
+  before_script:
+    - source .local/bin/env.sh
+    - ls -lrtha 
+    - ln -s "${LOCAL_DATA_PATH}/sim_output" sim_output
+    - ln -s "${LOCAL_DATA_PATH}/datasets/data" data
+    - ls -lrtha
+    - bash bin/get_calibrations
+
 include:
   - local: 'benchmarks/dis/config.yml'
-  - local: 'benchmarks/dvmp/config.yml'
+    #- local: 'benchmarks/dvmp/config.yml'
   - local: 'benchmarks/dvcs/config.yml'
 
 summary:
diff --git a/benchmarks/dvcs/config.yml b/benchmarks/dvcs/config.yml
index 3256516c10740dd96c6c33645d2967a2f8b58a08..63d78140a8fee8660eb340b3fc9114a036136d82 100644
--- a/benchmarks/dvcs/config.yml
+++ b/benchmarks/dvcs/config.yml
@@ -1,13 +1,17 @@
 dvcs:process:
   stage: process
-  timeout: 2 hour
+  extends: .phy_benchmark
+  tags:
+    - s3
   needs: ["common:detector"]
   script:
+    - bash benchmarks/dvcs/dvcs.sh --data-init
     - bash benchmarks/dvcs/dvcs.sh
 
 dvcs:results:
   stage: collect
   needs: ["dvcs:process"]
   script:
-    - pip install junitparser
+    - ls -lrth
+      #pip install junitparser
       #- python dvcs/scripts/merge_results.py
diff --git a/benchmarks/dvcs/dvcs.sh b/benchmarks/dvcs/dvcs.sh
index cc4b2a29a84e27c415e760f8c891ed6357f5fc90..bcbfa38cabc7c20c0ff89070183e4d0ebf9c15f0 100644
--- a/benchmarks/dvcs/dvcs.sh
+++ b/benchmarks/dvcs/dvcs.sh
@@ -1,28 +1,68 @@
 #!/bin/bash
 
+function print_the_help {
+  echo "USAGE: ${0} [--data-init]  "
+  echo "  OPTIONS: "
+  exit 
+}
+
+DATA_INIT=
+REC_ONLY=
+ANALYSIS_ONLY=
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    --data-init)
+      DATA_INIT=1
+      shift # past value
+      ;;
+    *)    # unknown option
+      #POSITIONAL+=("$1") # save it in an array for later
+      echo "unknown option $1"
+      print_the_help
+      shift # past argument
+      ;;
+  esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
 # these variables might not need exported.
-export JUGGLER_FILE_NAME_TAG="dvcs"
+export FILE_NAME_TAG="dvcs"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
-export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${FILE_NAME_TAG}.root"
 
 echo "JUGGLER_N_EVENTS      = ${JUGGLER_N_EVENTS}"
 echo "JUGGLER_DETECTOR      = ${JUGGLER_DETECTOR}"
-echo "JUGGLER_FILE_NAME_TAG = ${JUGGLER_FILE_NAME_TAG}"
+echo "FILE_NAME_TAG = ${FILE_NAME_TAG}"
+
+print_env.sh
 
 ## To run the reconstruction, we need the following global variables:
 ## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
 ## - JUGGLER_DETECTOR:       the detector package we want to use for this benchmark
 ## - DETECTOR_PATH:          full path to the detector definitions
 
-
-curl -o test_proton_dvcs_eic.hepmc "https://eicweb.phy.anl.gov/api/v4/projects/345/jobs/artifacts/master/raw/data/test_proton_dvcs_eic.hepmc?job=compile"
-if [[ "$?" -ne "0" ]] ; then
-  echo "Failed to download hepmc file"
-  exit 1
+if [[ -n "${DATA_INIT}" ]] ; then 
+  mc -C . config host add S3 https://dtn01.sdcc.bnl.gov:9000 $S3_ACCESS_KEY $S3_SECRET_KEY
+  mc -C . cat  --insecure S3/eictest/ATHENA/EVGEN/DVCS/DVCS_10x100_2M/DVCS.1.hepmc |  head  -n 1004 > "${LOCAL_DATA_PATH}/dvcs_test.hepmc"
+  if [[ "$?" -ne "0" ]] ; then
+    echo "Failed to download hepmc file"
+    exit 1
+  fi
+  exit
 fi
 
-export JUGGLER_N_EVENTS=10
+
+#curl -o test_proton_dvcs_eic.hepmc "https://eicweb.phy.anl.gov/api/v4/projects/345/jobs/artifacts/master/raw/data/test_proton_dvcs_eic.hepmc?job=compile"
+
 
 ## run geant4 simulations
 npsim --runType batch \
@@ -30,7 +70,7 @@ npsim --runType batch \
       -v ERROR \
       --numberOfEvents ${JUGGLER_N_EVENTS} \
       --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
-      --inputFiles test_proton_dvcs_eic.hepmc \
+      --inputFiles "${LOCAL_DATA_PATH}/dvcs_test.hepmc" \
       --outputFile  ${JUGGLER_SIM_FILE}
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running npsim"
@@ -46,6 +86,7 @@ if [[ "$?" -ne "0" ]] ; then
 fi
 
 mkdir -p results/dvcs
+rootls -t ${JUGGLER_SIM_FILE}
 
 root -b -q "benchmarks/dvcs/scripts/dvcs_tests.cxx(\"${JUGGLER_REC_FILE}\")"
 if [[ "$?" -ne "0" ]] ; then
diff --git a/benchmarks/dvcs/scripts/dvcs_ps_gen.cxx b/benchmarks/dvcs/scripts/dvcs_ps_gen.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c0642a7c05cd487a5668741ed00ca997d26c782f
--- /dev/null
+++ b/benchmarks/dvcs/scripts/dvcs_ps_gen.cxx
@@ -0,0 +1,33 @@
+void dvcs_ps_gen() {
+  double E_p = 100.0;
+  double M_p = 0.938;
+  double E_e = 5.0;
+   TLorentzVector target(0.0, 0.0, std::sqrt(E_p*E_p - M_p*M_p),E_p);
+   TLorentzVector beam(0.0, 0.0, E_e, std::sqrt(E_e*E_e+0.000511*0.000511));
+   TLorentzVector W = beam + target;
+
+   //(Momentum, Energy units are Gev/C, GeV)
+   Double_t masses[3] = { 0.000511,0.938, 0.0} ;
+
+   TGenPhaseSpace event;
+   if(!event.SetDecay(W, 3, masses))
+     std::cout << "derp\n";
+     ;
+
+   TH2F *h2 = new TH2F("h2","h2; Q^{2} ; t", 100,0,5, 100,-0.120,0);
+
+   for (Int_t n=0;n<1000000;n++) {
+      Double_t weight = event.Generate();
+
+      TLorentzVector* pElectron = event.GetDecay(0);
+      TLorentzVector* pProton = event.GetDecay(1);
+      TLorentzVector* pGamma = event.GetDecay(2);
+
+      TLorentzVector pq = beam - *pElectron;
+      TLorentzVector Delta = target - *pProton;
+
+      h2->Fill(-1.0*(pq.M2()) , Delta.M2() ,weight);
+      //std::cout <<  -1.0*(pq.M2())  << " , " <<  Delta.M2() << " , " << weight << "\n";
+   }
+   h2->Draw("colz");
+}
diff --git a/benchmarks/dvcs/scripts/dvcs_tests.cxx b/benchmarks/dvcs/scripts/dvcs_tests.cxx
index 8364dffc02d5c90ea04f63b53e7d2cc570eeef65..d7d0df8dade34e1b77d98ad8a72d43302392f1cd 100644
--- a/benchmarks/dvcs/scripts/dvcs_tests.cxx
+++ b/benchmarks/dvcs/scripts/dvcs_tests.cxx
@@ -106,7 +106,7 @@ void dvcs_tests(const char* fname = "rec_dvcs.root"){
   auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
                  .Define("thrownParticles", "mcparticles2[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
-                 .Define("dumRec", dumfourvec, {"DummyReconstructedParticles"})
+                 .Define("dumRec", dumfourvec, {"ReconstructedParticles"})
                  .Define("dumNPart", "dumRec.size()")
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("nTracks", "outputTrackParameters.size()")
diff --git a/benchmarks/dvcs/scripts/merge_results.py b/benchmarks/dvcs/scripts/merge_results.py
deleted file mode 100755
index 3ec0cca730710239bcb9d674934a0d0dc0850d50..0000000000000000000000000000000000000000
--- a/benchmarks/dvcs/scripts/merge_results.py
+++ /dev/null
@@ -1,64 +0,0 @@
-#!/usr/bin/env python3
-
-import json
-from junitparser import TestCase, TestSuite, JUnitXml, Skipped, Error, IntAttr, FloatAttr
-
-
-# Create the new element by subclassing Element or one of its child class,
-## and add custom attributes to it.
-#class MyTestCase(TestCase):
-#    foo = Attr()
-
-# Add the custom attribute
-#TestCase.id = IntAttr('id')
-TestCase.efficiency = FloatAttr('efficiency')
-#TestCase.custom = Attr('custom')
-#case = TestCase()
-#case.id = 123
-#case.rate = 0.95
-#case.custom = 'foobar'
-
-# After looking at two different python libraries (junit-xml and junitparser)
-# junitparser looks the most robust
-# https://github.com/weiwei/junitparser 
-
-def merge_results():
-    results = None;
-    with open("results/dvcs/dvcs_tests.json","r") as f: 
-        results = json.load(f)
-
-    # Create suite and add cases
-    suite = TestSuite('dvcs')
-    suite.add_property('energy', '10-on-100')
-
-    for tname,tres in results.items():
-        for ttype, tval in tres.items():
-            # Create cases
-            case1 = TestCase(tname)
-            case1.time = 1.0
-            case1.efficiency = tval
-            case1.classname = ttype
-            suite.add_testcase(case1)
-
-    xml = JUnitXml()
-    xml.add_testsuite(suite)
-    xml.write('results/dvcs/dvcs_report.xml',pretty=True)
-
-
-#test code for junit-xml:
-#from junit_xml import TestSuite, TestCase
-#    test_cases =  []
-#        print(test_name)
-#        print(test_res)
-#        for test_type, test_val in test_res.items():
-#           test_cases.append(TestCase(test_name, "dvcs.dvcs_tests.{}".format(test_type), 10, str(test_val), 'I am stderr!'))
-#    ts = TestSuite("my test suite", test_cases)
-#    # pretty printing is on by default but can be disabled using prettyprint=False
-#    print(TestSuite.to_xml_string([ts]))
-#    # you can also write the XML to a file and not pretty print it
-#    with open('results/dvcs/dvcs_report.xml', 'w') as f:
-#        TestSuite.to_file(f, [ts], prettyprint=True)
-
-if __name__ == "__main__":
-    # execute only if run as a script
-    merge_results()
diff --git a/benchmarks/dvmp/analysis/vm_invar.cxx b/benchmarks/dvmp/analysis/vm_invar.cxx
index fec278a549fcd40b55f2888ae69550646621d6e2..06fa8e15663c66f8bee7c7de9aff7cf1330d26f6 100644
--- a/benchmarks/dvmp/analysis/vm_invar.cxx
+++ b/benchmarks/dvmp/analysis/vm_invar.cxx
@@ -82,7 +82,7 @@ int vm_invar(const std::string& config_name)
   //====================================================================
 
   // Define analysis flow
-  auto d_im = d.Define("p_rec_sorted", momenta_sort_rec, {"DummyReconstructedParticles"})
+  auto d_im = d.Define("p_rec_sorted", momenta_sort_rec, {"ReconstructedParticles"})
                   .Define("p_sim_sorted", momenta_sort_sim, {"mcparticles2"})
                   .Define("N", "p_rec_sorted.size()")
                   .Define("invariant_quantities_rec", util::calc_inv_quant, {"p_rec_sorted"})
diff --git a/benchmarks/dvmp/analysis/vm_mass.cxx b/benchmarks/dvmp/analysis/vm_mass.cxx
index 9baeda4adbbe1905f86bfc1d2181a889bb03823c..0b90eaa61ce2651b97f3a6174172a346acc6ec73 100644
--- a/benchmarks/dvmp/analysis/vm_mass.cxx
+++ b/benchmarks/dvmp/analysis/vm_mass.cxx
@@ -96,7 +96,7 @@ int vm_mass(const std::string& config_name)
 
   // common_bench::PrintGeant4(mcparticles2);
   // Define analysis flow
-  auto d_im = d.Define("p_rec", common_bench::momenta_RC, {"DummyReconstructedParticles"})       //using dummy rc
+  auto d_im = d.Define("p_rec", common_bench::momenta_RC, {"ReconstructedParticles"})       //using dummy rc
                   .Define("N", "p_rec.size()")
                   .Define("p_sim", common_bench::momenta_from_simulation, {"mcparticles2"})
                   .Define("decay_pair_rec", find_decay_pair, {"p_rec"})
diff --git a/benchmarks/dvmp/config.yml b/benchmarks/dvmp/config.yml
index 32e929f9dd3e2648a047b1a9433ed9da98ff7cb5..57908c30671d00b469f10ade9dd342f863cd8d5d 100644
--- a/benchmarks/dvmp/config.yml
+++ b/benchmarks/dvmp/config.yml
@@ -1,6 +1,7 @@
 dvmp:generate:
   needs: ["common:detector"]
   image: eicweb.phy.anl.gov:4567/monte_carlo/lager/lager:unstable
+  extends: .phy_benchmark
   stage: generate
   timeout: 1 hours
   script:
@@ -12,6 +13,7 @@ dvmp:generate:
 
 dvmp:process:
   stage: process
+  extends: .phy_benchmark
   needs: ["common:detector", "dvmp:generate"]
   timeout: 2 hour
   script:
diff --git a/benchmarks/dvmp/dvmp.sh b/benchmarks/dvmp/dvmp.sh
index ac1c0af9693f9c5d02ad6828ddbe1148b563b02f..58a289fb130735a1c6931eaa2d4c9279d8a75985 100755
--- a/benchmarks/dvmp/dvmp.sh
+++ b/benchmarks/dvmp/dvmp.sh
@@ -85,8 +85,8 @@ echo "Running the digitization and reconstruction"
 export JUGGLER_SIM_FILE=${SIM_FILE}
 export JUGGLER_REC_FILE=${REC_FILE}
 xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
-  gaudirun.py options/reconstruction.py \
-  2>&1 > ${REC_LOG}
+  gaudirun.py options/reconstruction.py 
+##  2>&1 > ${REC_LOG}
 ## on-error, first retry running juggler again as there is still a random
 ## crash we need to address FIXME
 if [ "$?" -ne "0" ] ; then
diff --git a/bin/get_calibrations b/bin/get_calibrations
new file mode 100644
index 0000000000000000000000000000000000000000..57389c8e97117f9517ffde0a83823a9ac31dcc19
--- /dev/null
+++ b/bin/get_calibrations
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+branch=${1:-master}
+detector_benchmarks=https://eicweb.phy.anl.gov/EIC/benchmarks/detector_benchmarks/-/jobs/artifacts/${branch}/raw/
+
+mkdir -p config
+for i in results/emcal_barrel_calibration.json ; do
+  curl --fail -sL ${detector_benchmarks}/${i}?job=deploy_results --output config/$(basename ${i})
+done
diff --git a/options/reconstruction.py b/options/reconstruction.py
index 03b1ada35a4555d9b6fb4ee4a92af6a403b1868e..ca61cbaec805b5ebc55e06f6272cb2464b399598 100644
--- a/options/reconstruction.py
+++ b/options/reconstruction.py
@@ -4,160 +4,143 @@ from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
 from GaudiKernel import SystemOfUnits as units
 from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
 
+import json
 
-detector_name = "topside"
+detector_name = "athena"
 if "JUGGLER_DETECTOR" in os.environ :
   detector_name = str(os.environ["JUGGLER_DETECTOR"])
 
-# todo add checks
-input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
-output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
-n_events = str(os.environ["JUGGLER_N_EVENTS"])
-
-detector_path = detector_name
+detector_path = ""
 if "DETECTOR_PATH" in os.environ :
-    detector_path = str(os.environ["DETECTOR_PATH"])
+  detector_path = str(os.environ["DETECTOR_PATH"])
+
+compact_path = os.path.join(detector_path, detector_name)
 
-# get sampling fractions from system environment variable, 1.0 by default
+# RICH reconstruction
+qe_data = [(1.0, 0.25), (7.5, 0.25),]
+
+# CAL reconstruction
+# get sampling fractions from system environment variable
 ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253))
-cb_ecal_sf = float(os.environ.get("CB_ECAL_SAMP_FRAC", 0.01324))
 cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 0.038))
 ci_hcal_sf = float(os.environ.get("CI_HCAL_SAMP_FRAC", 0.025))
 ce_hcal_sf = float(os.environ.get("CE_HCAL_SAMP_FRAC", 0.025))
-scifi_barrel_sf = float(os.environ.get("CB_EMCAL_SCFI_SAMP_FRAC",0.0938))
 
+# input arguments from calibration file
+with open('config/emcal_barrel_calibration.json') as f:
+    calib_data = json.load(f)['electron']
+
+print(calib_data)
 
-geo_service  = GeoSvc("GeoSvc",
-        detectors=["{}/{}.xml".format(detector_path, detector_name)])
-podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=DEBUG)
+img_barrel_sf = float(calib_data['sampling_fraction_img'])
+scifi_barrel_sf = float(calib_data['sampling_fraction_scfi'])
 
+# input and output
+input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
+output_rec = str(os.environ["JUGGLER_REC_FILE"])
+n_events = int(os.environ["JUGGLER_N_EVENTS"])
+
+# geometry service
+geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)], OutputLevel=WARNING)
+# data service
+podioevent = EICDataSvc("EventDataSvc", inputs=input_sims, OutputLevel=WARNING)
+
+# juggler components
 from Configurables import PodioInput
 from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
 from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
 from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__PhotoMultiplierHitCollection_dd4pod__PhotoMultiplierHitCollection_ as PMTCopier
 
-from Configurables import Jug__Digi__SiliconTrackerDigi as TrackerDigi
-
-from Configurables import Jug__Base__MC2DummyParticle as MC2DummyParticle
+from Configurables import Jug__Digi__PhotoMultiplierDigi as PhotoMultiplierDigi
+from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
+from Configurables import Jug__Digi__UFSDTrackerDigi as TrackerDigi
 
 from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
-
+from Configurables import Jug__Reco__TrackingHitsCollector2 as TrackingHitsCollector
 from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
-from Configurables import Jug__Reco__Tracker2SourceLinker as Tracker2SourceLinker
-#from Configurables import Jug__Reco__TrackerSourcesLinker as TrackerSourcesLinker
+from Configurables import Jug__Reco__TrackerSourcesLinker as TrackerSourcesLinker
 #from Configurables import Jug__Reco__TrackingHitsSourceLinker as TrackingHitsSourceLinker
 from Configurables import Jug__Reco__TrackParamTruthInit as TrackParamTruthInit
 from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
 from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
-
 from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
 from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
 
-from Configurables import Jug__Digi__EMCalorimeterDigi as EMCalorimeterDigi
-from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
-from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
 from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
-from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
-from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
 from Configurables import Jug__Reco__CalorimeterHitsMerger as CalHitsMerger
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+
 from Configurables import Jug__Reco__ImagingPixelReco as ImCalPixelReco
 from Configurables import Jug__Reco__ImagingTopoCluster as ImagingCluster
-from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
-
-
-
-
 
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
 
-podioinput = PodioInput("PodioReader",
-        collections=["mcparticles",
-            "TrackerEndcapHits","TrackerBarrelHits",
-            "VertexBarrelHits","VertexEndcapHits",
-            "EcalEndcapNHits", "EcalEndcapPHits",
-            "EcalBarrelHits", "EcalBarrelScFiHits", 
-            "HcalEndcapPHits", "HcalEndcapNHits",
-            "HcalBarrelHits",
-            ])#, OutputLevel=DEBUG)
+from Configurables import Jug__Reco__PhotoMultiplierReco as PhotoMultiplierReco
+from Configurables import Jug__Reco__PhotoRingClusters as PhotoRingClusters
 
-dummy = MC2DummyParticle("MC2Dummy",
-        inputCollection="mcparticles",
-        outputCollection="DummyReconstructedParticles",
-        smearing = 0.1)
+from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
 
+# branches needed from simulation root file
+sim_coll = [
+    "mcparticles",
+    "EcalEndcapNHits",
+    "EcalEndcapPHits",
+    "EcalBarrelHits",
+    "EcalBarrelScFiHits",
+    "HcalBarrelHits",
+    "HcalEndcapPHits",
+    "HcalEndcapNHits",
+    "TrackerEndcapHits",
+    "TrackerBarrelHits",
+    "GEMTrackerEndcapHits",
+    "VertexBarrelHits",
+    "VertexEndcapHits",
+    "DRICHHits",
+    "MRICHHits"
+]
+
+# list of algorithms
+algorithms = []
+
+# input
+podin = PodioInput("PodioReader", collections=sim_coll)
+algorithms.append(podin)
 
 ## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
-copier = MCCopier("MCCopier", 
-        inputCollection="mcparticles", 
-        outputCollection="mcparticles2") 
-trkcopier = TrkCopier("TrkCopier", 
-        inputCollection="TrackerBarrelHits", 
-        outputCollection="TrackerBarrelHits2") 
-algorithms = [podioinput,dummy, copier, trkcopier]                                                                              
-
-# Tracker and vertex digitization
-trk_b_digi = TrackerDigi("trk_b_digi", 
-        inputHitCollection="TrackerBarrelHits",
-        outputHitCollection="TrackerBarrelRawHits",
-        timeResolution=8)
-algorithms.append(trk_b_digi)
-
-trk_ec_digi = TrackerDigi("trk_ec_digi", 
-        inputHitCollection="TrackerEndcapHits",
-        outputHitCollection="TrackerEndcapRawHits",
-        timeResolution=8)
-algorithms.append(trk_ec_digi)
-
-vtx_b_digi = TrackerDigi("vtx_b_digi", 
-        inputHitCollection="VertexBarrelHits",
-        outputHitCollection="VertexBarrelRawHits",
-        timeResolution=8)
-algorithms.append(vtx_b_digi)
-
-vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
-        inputHitCollection="VertexEndcapHits",
-        outputHitCollection="VertexEndcapRawHits",
-        timeResolution=8)
-algorithms.append(vtx_ec_digi)
-
-# Tracker and vertex reconstruction
-trk_b_reco = TrackerHitReconstruction("trk_b_reco",
-        inputHitCollection = trk_b_digi.outputHitCollection,
-        outputHitCollection="TrackerBarrelRecHits")
-algorithms.append(trk_b_reco)
-
-trk_ec_reco = TrackerHitReconstruction("trk_ec_reco",
-        inputHitCollection = trk_ec_digi.outputHitCollection,
-        outputHitCollection="TrackerEndcapRecHits")
-algorithms.append(trk_ec_reco)
-
-vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
-        inputHitCollection = vtx_b_digi.outputHitCollection,
-        outputHitCollection="VertexBarrelRecHits")
-algorithms.append(vtx_b_reco)
+copier = MCCopier("MCCopier",
+        inputCollection="mcparticles",
+        outputCollection="mcparticles2")
+algorithms.append(copier)
 
-vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
-        inputHitCollection = vtx_ec_digi.outputHitCollection,
-        outputHitCollection="VertexEndcapRecHits")
-algorithms.append(vtx_ec_reco)
+trkcopier = TrkCopier("TrkCopier",
+        inputCollection="TrackerBarrelHits",
+        outputCollection="TrackerBarrelHits2")
+algorithms.append(trkcopier)
 
+pmtcopier = PMTCopier("PMTCopier",
+        inputCollection="DRICHHits",
+        outputCollection="DRICHHits2")
+algorithms.append(pmtcopier)
 
 # Crystal Endcap Ecal
 ce_ecal_daq = dict(
-        dynamicRangeADC=5.*GeV,
+        dynamicRangeADC=5.*units.GeV,
         capacityADC=32768,
         pedestalMean=400,
         pedestalSigma=3)
 
 ce_ecal_digi = CalHitDigi("ce_ecal_digi",
         inputHitCollection="EcalEndcapNHits",
-        outputHitCollection="EcalEndcapNHitsDigi",
+        outputHitCollection="EcalEndcapNRawHits",
         energyResolutions=[0., 0.02, 0.],
         **ce_ecal_daq)
 algorithms.append(ce_ecal_digi)
 
 ce_ecal_reco = CalHitReco("ce_ecal_reco",
         inputHitCollection=ce_ecal_digi.outputHitCollection,
-        outputHitCollection="EcalEndcapNHitsReco",
+        outputHitCollection="EcalEndcapNRecHits",
         thresholdFactor=4,          # 4 sigma cut on pedestal sigma
         readoutClass="EcalEndcapNHits",
         sectorField="sector",
@@ -165,62 +148,58 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco",
 algorithms.append(ce_ecal_reco)
 
 ce_ecal_cl = IslandCluster("ce_ecal_cl",
-        # OutputLevel=DEBUG,
         inputHitCollection=ce_ecal_reco.outputHitCollection,
         outputProtoClusterCollection="EcalEndcapNProtoClusters",
         splitCluster=False,
-        minClusterHitEdep=1.0*MeV,  # discard low energy hits
-        minClusterCenterEdep=30*MeV,
-        sectorDist=5.0*cm,
-        dimScaledLocalDistXY=[1.8, 1.8])          # dimension scaled dist is good for hybrid sectors with different module size
+        minClusterHitEdep=1.0*units.MeV,  # discard low energy hits
+        minClusterCenterEdep=30*units.MeV,
+        sectorDist=5.0*units.cm,
+        dimScaledLocalDistXY=[1.8, 1.8]) # dimension scaled dist is good for hybrid sectors with different module size
 algorithms.append(ce_ecal_cl)
 
 ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
         inputHitCollection=ce_ecal_cl.inputHitCollection,
         inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
         outputClusterCollection="EcalEndcapNClusters",
-        outputInfoCollection="EcalEndcapNClustersInfo",
         samplingFraction=0.998,      # this accounts for a small fraction of leakage
         logWeightBase=4.6)
 algorithms.append(ce_ecal_clreco)
 
 # Endcap Sampling Ecal
 ci_ecal_daq = dict(
-        dynamicRangeADC=50.*MeV,
+        dynamicRangeADC=50.*units.MeV,
         capacityADC=32768,
         pedestalMean=400,
         pedestalSigma=10)
 
 ci_ecal_digi = CalHitDigi("ci_ecal_digi",
         inputHitCollection="EcalEndcapPHits",
-        outputHitCollection="EcalEndcapPHitsDigi",
+        outputHitCollection="EcalEndcapPRawHits",
         **ci_ecal_daq)
 algorithms.append(ci_ecal_digi)
 
 ci_ecal_reco = CalHitReco("ci_ecal_reco",
         inputHitCollection=ci_ecal_digi.outputHitCollection,
-        outputHitCollection="EcalEndcapPHitsReco",
+        outputHitCollection="EcalEndcapPRecHits",
         thresholdFactor=5.0,
         **ci_ecal_daq)
 algorithms.append(ci_ecal_reco)
 
 # merge hits in different layer (projection to local x-y plane)
 ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
-        # OutputLevel=DEBUG,
         inputHitCollection=ci_ecal_reco.outputHitCollection,
-        outputHitCollection="EcalEndcapPHitsRecoXY",
+        outputHitCollection="EcalEndcapPRecHitsXY",
         fields=["layer", "slice"],
         fieldRefNumbers=[1, 0],
         readoutClass="EcalEndcapPHits")
 algorithms.append(ci_ecal_merger)
 
 ci_ecal_cl = IslandCluster("ci_ecal_cl",
-        # OutputLevel=DEBUG,
         inputHitCollection=ci_ecal_merger.outputHitCollection,
-        outputProtoClusterCollection="EcalEndcapProtoClusters",
+        outputProtoClusterCollection="EcalEndcapPProtoClusters",
         splitCluster=False,
-        minClusterCenterEdep=10.*MeV,
-        localDistXY=[10*mm, 10*mm])
+        minClusterCenterEdep=10.*units.MeV,
+        localDistXY=[10*units.mm, 10*units.mm])
 algorithms.append(ci_ecal_cl)
 
 ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
@@ -233,49 +212,48 @@ ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
 algorithms.append(ci_ecal_clreco)
 
 # Central Barrel Ecal (Imaging Cal.)
-cb_ecal_daq = dict(
-        dynamicRangeADC=3*MeV,
+img_barrel_daq = dict(
+        dynamicRangeADC=3*units.MeV,
         capacityADC=8192,
         pedestalMean=400,
         pedestalSigma=20)   # about 6 keV
 
-cb_ecal_digi = CalHitDigi("cb_ecal_digi",
+img_barrel_digi = CalHitDigi("img_barrel_digi",
         inputHitCollection="EcalBarrelHits",
-        outputHitCollection="EcalBarrelHitsDigi",
+        outputHitCollection="EcalBarrelImagingRawHits",
         energyResolutions=[0., 0.02, 0.],   # 2% flat resolution
-        **cb_ecal_daq)
-algorithms.append(cb_ecal_digi)
+        **img_barrel_daq)
+algorithms.append(img_barrel_digi)
 
-cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
-        inputHitCollection=cb_ecal_digi.outputHitCollection,
-        outputHitCollection="EcalBarrelHitsReco",
+img_barrel_reco = ImCalPixelReco("img_barrel_reco",
+        inputHitCollection=img_barrel_digi.outputHitCollection,
+        outputHitCollection="EcalBarrelImagingRecHits",
         thresholdFactor=3,  # about 20 keV
         readoutClass="EcalBarrelHits",  # readout class
         layerField="layer",             # field to get layer id
         sectorField="module",           # field to get sector id
-        **cb_ecal_daq)
-algorithms.append(cb_ecal_reco)
-
-cb_ecal_cl = ImagingCluster("cb_ecal_cl",
-        inputHitCollection=cb_ecal_reco.outputHitCollection,
-        outputProtoClusterCollection="EcalBarrelProtoClusters",
-        localDistXY=[2.*mm, 2*mm],              # same layer
-        layerDistEtaPhi=[10*mrad, 10*mrad],     # adjacent layer
+        **img_barrel_daq)
+algorithms.append(img_barrel_reco)
+
+img_barrel_cl = ImagingCluster("img_barrel_cl",
+        inputHitCollection=img_barrel_reco.outputHitCollection,
+        outputProtoClusterCollection="EcalBarrelImagingProtoClusters",
+        localDistXY=[2.*units.mm, 2*units.mm],  # same layer
+        layerDistEtaPhi=[10*units.mrad, 10*units.mrad],     # adjacent layer
         neighbourLayersRange=2,                 # id diff for adjacent layer
-        sectorDist=3.*cm)                       # different sector
-algorithms.append(cb_ecal_cl)
-
-cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
-        samplingFraction=cb_ecal_sf,
-        inputHitCollection=cb_ecal_cl.inputHitCollection,
-        inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection,
-        outputClusterCollection="EcalBarrelClusters",
-        outputInfoCollection="EcalBarrelClustersInfo",
-        outputLayerCollection="EcalBarrelLayers")
-algorithms.append(cb_ecal_clreco)
-
-#Central ECAL SciFi
-# use the same daq_setting for digi/reco pair
+        sectorDist=3.*units.cm)                       # different sector
+algorithms.append(img_barrel_cl)
+
+img_barrel_clreco = ImagingClusterReco("img_barrel_clreco",
+        samplingFraction=img_barrel_sf,
+        inputHitCollection=img_barrel_cl.inputHitCollection,
+        inputProtoClusterCollection=img_barrel_cl.outputProtoClusterCollection,
+        outputClusterCollection="EcalBarrelImagingClusters",
+        outputInfoCollection="EcalBarrelImagingClustersInfo",
+        outputLayerCollection="EcalBarrelImagingLayers")
+algorithms.append(img_barrel_clreco)
+
+# Central ECAL SciFi
 scfi_barrel_daq = dict(
         dynamicRangeADC=50.*MeV,
         capacityADC=32768,
@@ -284,13 +262,13 @@ scfi_barrel_daq = dict(
 
 scfi_barrel_digi = CalHitDigi("scfi_barrel_digi",
         inputHitCollection="EcalBarrelScFiHits",
-        outputHitCollection="EcalBarrelScFiHitsDigi",
+        outputHitCollection="EcalBarrelScFiRawHits",
         **scfi_barrel_daq)
 algorithms.append(scfi_barrel_digi)
 
 scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
         inputHitCollection=scfi_barrel_digi.outputHitCollection,
-        outputHitCollection="EcalBarrelScFiHitsReco",
+        outputHitCollection="EcalBarrelScFiRecHits",
         thresholdFactor=5.0,
         readoutClass="EcalBarrelScFiHits",
         layerField="layer",
@@ -301,49 +279,46 @@ algorithms.append(scfi_barrel_reco)
 
 # merge hits in different layer (projection to local x-y plane)
 scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
-        # OutputLevel=DEBUG,
-        inputHitCollection=scfi_barrel_reco.outputHitCollection,
-        outputHitCollection="EcalBarrelScFiGridReco",
-        fields=["fiber"],
-        fieldRefNumbers=[1],
-        readoutClass="EcalBarrelScFiHits")
+         inputHitCollection=scfi_barrel_reco.outputHitCollection,
+         outputHitCollection="EcalBarrelScFiGridReco",
+         fields=["fiber"],
+         fieldRefNumbers=[1],
+         readoutClass="EcalBarrelScFiHits")
 algorithms.append(scfi_barrel_merger)
 
 scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
-        # OutputLevel=DEBUG,
-        inputHitCollection=scfi_barrel_merger.outputHitCollection,
-        outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
-        splitCluster=False,
-        minClusterCenterEdep=10.*MeV,
-        localDistXZ=[30*mm, 30*mm])
+         inputHitCollection=scfi_barrel_merger.outputHitCollection,
+         outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
+         splitCluster=False,
+         minClusterCenterEdep=10.*MeV,
+         localDistXZ=[30*mm, 30*mm])
 algorithms.append(scfi_barrel_cl)
 
 scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
-       inputHitCollection=scfi_barrel_cl.inputHitCollection,
-       inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection,
-       outputClusterCollection="EcalBarrelScFiClusters",
-       outputInfoCollection="EcalBarrelScFiClustersInfo",
-       logWeightBase=6.2,
-       samplingFraction= scifi_barrel_sf)
+         inputHitCollection=scfi_barrel_cl.inputHitCollection,
+         inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection,
+         outputClusterCollection="EcalBarrelScFiClusters",
+         outputInfoCollection="EcalBarrelScFiClustersInfo",
+         logWeightBase=6.2,
+         samplingFraction= scifi_barrel_sf)
 algorithms.append(scfi_barrel_clreco)
 
-
 # Central Barrel Hcal
 cb_hcal_daq = dict(
-         dynamicRangeADC=50.*MeV,
+         dynamicRangeADC=50.*units.MeV,
          capacityADC=32768,
          pedestalMean=400,
          pedestalSigma=10)
 
 cb_hcal_digi = CalHitDigi("cb_hcal_digi",
          inputHitCollection="HcalBarrelHits",
-         outputHitCollection="HcalBarrelHitsDigi",
+         outputHitCollection="HcalBarrelRawHits",
          **cb_hcal_daq)
 algorithms.append(cb_hcal_digi)
 
 cb_hcal_reco = CalHitReco("cb_hcal_reco",
         inputHitCollection=cb_hcal_digi.outputHitCollection,
-        outputHitCollection="HcalBarrelHitsReco",
+        outputHitCollection="HcalBarrelRecHits",
         thresholdFactor=5.0,
         readoutClass="HcalBarrelHits",
         layerField="layer",
@@ -353,7 +328,7 @@ algorithms.append(cb_hcal_reco)
 
 cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
         inputHitCollection=cb_hcal_reco.outputHitCollection,
-        outputHitCollection="HcalBarrelHitsRecoXY",
+        outputHitCollection="HcalBarrelRecHitsXY",
         readoutClass="HcalBarrelHits",
         fields=["layer", "slice"],
         fieldRefNumbers=[1, 0])
@@ -363,8 +338,8 @@ cb_hcal_cl = IslandCluster("cb_hcal_cl",
         inputHitCollection=cb_hcal_merger.outputHitCollection,
         outputProtoClusterCollection="HcalBarrelProtoClusters",
         splitCluster=False,
-        minClusterCenterEdep=30.*MeV,
-        localDistXY=[15.*cm, 15.*cm])
+        minClusterCenterEdep=30.*units.MeV,
+        localDistXY=[15.*units.cm, 15.*units.cm])
 algorithms.append(cb_hcal_cl)
 
 cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
@@ -378,26 +353,27 @@ algorithms.append(cb_hcal_clreco)
 
 # Hcal Hadron Endcap
 ci_hcal_daq = dict(
-         dynamicRangeADC=50.*MeV,
+         dynamicRangeADC=50.*units.MeV,
          capacityADC=32768,
          pedestalMean=400,
          pedestalSigma=10)
+
 ci_hcal_digi = CalHitDigi("ci_hcal_digi",
          inputHitCollection="HcalEndcapPHits",
-         outputHitCollection="HcalEndcapPHitsDigi",
+         outputHitCollection="HcalEndcapPRawHits",
          **ci_hcal_daq)
 algorithms.append(ci_hcal_digi)
 
 ci_hcal_reco = CalHitReco("ci_hcal_reco",
         inputHitCollection=ci_hcal_digi.outputHitCollection,
-        outputHitCollection="HcalEndcapPHitsReco",
+        outputHitCollection="HcalEndcapPRecHits",
         thresholdFactor=5.0,
         **ci_hcal_daq)
 algorithms.append(ci_hcal_reco)
 
 ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
         inputHitCollection=ci_hcal_reco.outputHitCollection,
-        outputHitCollection="HcalEndcapPHitsRecoXY",
+        outputHitCollection="HcalEndcapPRecHitsXY",
         readoutClass="HcalEndcapPHits",
         fields=["layer", "slice"],
         fieldRefNumbers=[1, 0])
@@ -407,8 +383,8 @@ ci_hcal_cl = IslandCluster("ci_hcal_cl",
         inputHitCollection=ci_hcal_merger.outputHitCollection,
         outputProtoClusterCollection="HcalEndcapPProtoClusters",
         splitCluster=False,
-        minClusterCenterEdep=30.*MeV,
-        localDistXY=[15.*cm, 15.*cm])
+        minClusterCenterEdep=30.*units.MeV,
+        localDistXY=[15.*units.cm, 15.*units.cm])
 algorithms.append(ci_hcal_cl)
 
 ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
@@ -422,27 +398,27 @@ algorithms.append(ci_hcal_clreco)
 
 # Hcal Electron Endcap
 ce_hcal_daq = dict(
-        dynamicRangeADC=50.*MeV,
+        dynamicRangeADC=50.*units.MeV,
         capacityADC=32768,
         pedestalMean=400,
         pedestalSigma=10)
 
 ce_hcal_digi = CalHitDigi("ce_hcal_digi",
         inputHitCollection="HcalEndcapNHits",
-        outputHitCollection="HcalEndcapNHitsDigi",
+        outputHitCollection="HcalEndcapNRawHits",
         **ce_hcal_daq)
 algorithms.append(ce_hcal_digi)
 
 ce_hcal_reco = CalHitReco("ce_hcal_reco",
         inputHitCollection=ce_hcal_digi.outputHitCollection,
-        outputHitCollection="HcalEndcapNHitsReco",
+        outputHitCollection="HcalEndcapNRecHits",
         thresholdFactor=5.0,
         **ce_hcal_daq)
 algorithms.append(ce_hcal_reco)
 
 ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
         inputHitCollection=ce_hcal_reco.outputHitCollection,
-        outputHitCollection="HcalEndcapNHitsRecoXY",
+        outputHitCollection="HcalEndcapNRecHitsXY",
         readoutClass="HcalEndcapNHits",
         fields=["layer", "slice"],
         fieldRefNumbers=[1, 0])
@@ -452,8 +428,8 @@ ce_hcal_cl = IslandCluster("ce_hcal_cl",
         inputHitCollection=ce_hcal_merger.outputHitCollection,
         outputProtoClusterCollection="HcalEndcapNProtoClusters",
         splitCluster=False,
-        minClusterCenterEdep=30.*MeV,
-        localDistXY=[15.*cm, 15.*cm])
+        minClusterCenterEdep=30.*units.MeV,
+        localDistXY=[15.*units.cm, 15.*units.cm])
 algorithms.append(ce_hcal_cl)
 
 ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
@@ -465,119 +441,154 @@ ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
         samplingFraction=ce_hcal_sf)
 algorithms.append(ce_hcal_clreco)
 
-# Track Source linker 
-sourcelinker = TrackerSourceLinker("trk_srclinker",
-        inputHitCollection="TrackerBarrelRecHits",
-        outputSourceLinks="BarrelTrackSourceLinks",
-        OutputLevel=DEBUG)
-algorithms.append(sourcelinker)
+# Tracking
+trk_b_digi = TrackerDigi("trk_b_digi",
+        inputHitCollection="TrackerBarrelHits",
+        outputHitCollection="TrackerBarrelRawHits",
+        timeResolution=8)
+algorithms.append(trk_b_digi)
+
+trk_ec_digi = TrackerDigi("trk_ec_digi",
+        inputHitCollection="TrackerEndcapHits",
+        outputHitCollection="TrackerEndcapRawHits",
+        timeResolution=8)
+algorithms.append(trk_ec_digi)
+
+vtx_b_digi = TrackerDigi("vtx_b_digi",
+        inputHitCollection="VertexBarrelHits",
+        outputHitCollection="VertexBarrelRawHits",
+        timeResolution=8)
+algorithms.append(vtx_b_digi)
+
+vtx_ec_digi = TrackerDigi("vtx_ec_digi",
+        inputHitCollection="VertexEndcapHits",
+        outputHitCollection="VertexEndcapRawHits",
+        timeResolution=8)
+algorithms.append(vtx_ec_digi)
 
-trk_hits_srclnkr = Tracker2SourceLinker("trk_hits_srclnkr",
-        TrackerBarrelHits="TrackerBarrelRecHits",
-        TrackerEndcapHits="TrackerEndcapRecHits",
-        outputMeasurements="lnker2Measurements",
-        outputSourceLinks="lnker2Links",
-        allTrackerHits="linker2AllHits",
-        OutputLevel=DEBUG)
-algorithms.append(trk_hits_srclnkr)
+gem_ec_digi = TrackerDigi("gem_ec_digi",
+        inputHitCollection="GEMTrackerEndcapHits",
+        outputHitCollection="GEMTrackerEndcapRawHits",
+        timeResolution=10)
+algorithms.append(gem_ec_digi)
+
+# Tracker and vertex reconstruction
+trk_b_reco = TrackerHitReconstruction("trk_b_reco",
+        inputHitCollection = trk_b_digi.outputHitCollection,
+        outputHitCollection="TrackerBarrelRecHits")
+algorithms.append(trk_b_reco)
+
+trk_ec_reco = TrackerHitReconstruction("trk_ec_reco",
+        inputHitCollection = trk_ec_digi.outputHitCollection,
+        outputHitCollection="TrackerEndcapRecHits")
+algorithms.append(trk_ec_reco)
+
+vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
+        inputHitCollection = vtx_b_digi.outputHitCollection,
+        outputHitCollection="VertexBarrelRecHits")
+algorithms.append(vtx_b_reco)
+
+vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+        inputHitCollection = vtx_ec_digi.outputHitCollection,
+        outputHitCollection="VertexEndcapRecHits")
+algorithms.append(vtx_ec_reco)
+
+gem_ec_reco = TrackerHitReconstruction("gem_ec_reco",
+        inputHitCollection=gem_ec_digi.outputHitCollection,
+        outputHitCollection="GEMTrackerEndcapRecHits")
+algorithms.append(gem_ec_reco)
+
+# Tracking hit collector
+trk_hit_col = TrackingHitsCollector("trk_hit_col",
+        inputTrackingHits=[
+            str(trk_b_reco.outputHitCollection),
+            str(trk_ec_reco.outputHitCollection),
+            str(vtx_b_reco.outputHitCollection),
+            str(vtx_ec_reco.outputHitCollection),
+            str(gem_ec_reco.outputHitCollection) ],
+        trackingHits="trackingHits")
+algorithms.append(trk_hit_col)
+
+# Hit Source linker
+sourcelinker = TrackerSourceLinker("trk_srcslnkr",
+        inputHitCollection = trk_hit_col.trackingHits,
+        outputSourceLinks = "TrackSourceLinks",
+        outputMeasurements = "TrackMeasurements")
+algorithms.append(sourcelinker)
 
 ## Track param init
 truth_trk_init = TrackParamTruthInit("truth_trk_init",
         inputMCParticles="mcparticles",
-        outputInitialTrackParameters="InitTrackParams",
-        OutputLevel=DEBUG)
+        outputInitialTrackParameters="InitTrackParams")
 algorithms.append(truth_trk_init)
 
-#clust_trk_init = TrackParamClusterInit("clust_trk_init",
-#        inputClusters="SimpleClusters",
-#        outputInitialTrackParameters="InitTrackParamsFromClusters",
-#        OutputLevel=DEBUG)
-#algorithms.append(clust_trk_init)
-
-#vtxcluster_trk_init = TrackParamVertexClusterInit("vtxcluster_trk_init",
-#        inputVertexHits="VertexBarrelRecHits",
-#        inputClusters="SimpleClusters",
-#        outputInitialTrackParameters="InitTrackParamsFromVtxClusters",
-#        maxHitRadius=40.0*units.mm,
-#        OutputLevel=DEBUG)
-
 # Tracking algorithms
 trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
         inputSourceLinks = sourcelinker.outputSourceLinks,
         inputMeasurements = sourcelinker.outputMeasurements,
-        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
-        outputTrajectories="trajectories",
-        OutputLevel=DEBUG)
+        inputInitialTrackParameters = truth_trk_init.outputInitialTrackParameters,
+        outputTrajectories = "trajectories")
 algorithms.append(trk_find_alg)
 
 parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
-        inputTrajectories="trajectories",
-        outputParticles="ReconstructedParticles",
-        outputTrackParameters="outputTrackParameters",
-        OutputLevel=DEBUG)
+        inputTrajectories = trk_find_alg.outputTrajectories,
+        outputParticles = "ReconstructedParticles",
+        outputTrackParameters = "outputTrackParameters")
 algorithms.append(parts_from_fit)
 
-#trk_find_alg1 = TrackFindingAlgorithm("trk_find_alg1",
-#        inputSourceLinks = trk_hits_srclnkr.outputSourceLinks,
-#        inputMeasurements = trk_hits_srclnkr.outputMeasurements,
-#        inputInitialTrackParameters= "InitTrackParamsFromClusters", 
-#        outputTrajectories="trajectories1",
-#        OutputLevel=DEBUG)
-#parts_from_fit1 = ParticlesFromTrackFit("parts_from_fit1",
-#        inputTrajectories="trajectories1",
-#        outputParticles="ReconstructedParticles1",
-#        outputTrackParameters="outputTrackParameters1",
-#        OutputLevel=DEBUG)
-#
-#trk_find_alg2 = TrackFindingAlgorithm("trk_find_alg2",
-#        inputSourceLinks = trk_hits_srclnkr.outputSourceLinks,
-#        inputMeasurements = trk_hits_srclnkr.outputMeasurements,
-#        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
-#        #inputInitialTrackParameters= "InitTrackParamsFromVtxClusters", 
-#        outputTrajectories="trajectories2",
-#        OutputLevel=DEBUG)
-#parts_from_fit2 = ParticlesFromTrackFit("parts_from_fit2",
-#        inputTrajectories="trajectories2",
-#        outputParticles="ReconstructedParticles2",
-#        outputTrackParameters="outputTrackParameters2",
-#        OutputLevel=DEBUG)
-
-
-#types = []
-## this printout is useful to check that the type information is passed to python correctly
-#print("---------------------------------------\n")
-#print("---\n# List of input and output types by class")
-#for configurable in sorted([ PodioInput, EICDataSvc, PodioOutput,
-#                             TrackerHitReconstruction,ExampleCaloDigi, 
-#                             UFSDTrackerDigi, TrackerSourceLinker,
-#                             PodioOutput],
-#                           key=lambda c: c.getType()):
-#    print("\"{}\":".format(configurable.getType()))
-#    props = configurable.getDefaultProperties()
-#    for propname, prop in sorted(props.items()):
-#        print(" prop name: {}".format(propname))
-#        if isinstance(prop, DataHandleBase):
-#            types.append(prop.type())
-#            print("  {}: \"{}\"".format(propname, prop.type()))
-#print("---")
-
-out = PodioOutput("out", filename=output_rec_file)
-out.outputCommands = ["keep *", 
-        "drop BarrelTrackSourceLinks", 
+# DRICH
+drich_digi = PhotoMultiplierDigi("drich_digi",
+        inputHitCollection="DRICHHits",
+        outputHitCollection="DRICHRawHits",
+        quantumEfficiency=[(a*units.eV, b) for a, b in qe_data])
+algorithms.append(drich_digi)
+
+drich_reco = PhotoMultiplierReco("drich_reco",
+        inputHitCollection=drich_digi.outputHitCollection,
+        outputHitCollection="DRICHRecHits")
+algorithms.append(drich_reco)
+
+# FIXME
+#drich_cluster = PhotoRingClusters("drich_cluster",
+#        inputHitCollection=pmtreco.outputHitCollection,
+#        #inputTrackCollection="ReconstructedParticles",
+#        outputClusterCollection="ForwardRICHClusters")
+
+# MRICH
+mrich_digi = PhotoMultiplierDigi("mrich_digi",
+        inputHitCollection="MRICHHits",
+        outputHitCollection="MRICHRawHits",
+        quantumEfficiency=[(a*units.eV, b) for a, b in qe_data])
+algorithms.append(mrich_digi)
+
+mrich_reco = PhotoMultiplierReco("mrich_reco",
+        inputHitCollection=mrich_digi.outputHitCollection,
+        outputHitCollection="MRICHRecHits")
+algorithms.append(mrich_reco)
+
+# FIXME
+#mrich_cluster = PhotoRingClusters("drich_cluster",
+#        inputHitCollection=pmtreco.outputHitCollection,
+#        #inputTrackCollection="ReconstructedParticles",
+#        outputClusterCollection="ForwardRICHClusters")
+
+# Output
+podout = PodioOutput("out", filename=output_rec)
+podout.outputCommands = [
+        "keep *",
+        "drop *Digi",
+        "keep *Reco*",
+        "keep *ClusterHits",
+        "keep *Clusters",
+        "keep *Layers",
         "drop InitTrackParams",
-        "drop trajectories",
-        "drop outputSourceLinks",
-        "drop outputInitialTrackParameters",
-        "drop mcparticles"
-        ]
-algorithms.append(out)
+        ] + [ "drop " + c for c in sim_coll]
+algorithms.append(podout)
 
 ApplicationMgr(
     TopAlg = algorithms,
     EvtSel = 'NONE',
     EvtMax   = n_events,
     ExtSvc = [podioevent,geo_service],
-    OutputLevel=DEBUG
+    OutputLevel=WARNING
  )
-
-
diff --git a/options/tracker_reconstruction.py b/options/tracker_reconstruction.py
deleted file mode 100644
index 92a2a07c7403e732ef39065958c617376b862f96..0000000000000000000000000000000000000000
--- a/options/tracker_reconstruction.py
+++ /dev/null
@@ -1,238 +0,0 @@
-from Gaudi.Configuration import *
-
-from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
-from GaudiKernel import SystemOfUnits as units
-
-detector_name = "topside"
-if "JUGGLER_DETECTOR" in os.environ :
-  detector_name = str(os.environ["JUGGLER_DETECTOR"])
-
-# todo add checks
-input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
-output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
-n_events = str(os.environ["JUGGLER_N_EVENTS"])
-
-detector_path = detector_name
-if "DETECTOR_PATH" in os.environ :
-    detector_path = str(os.environ["DETECTOR_PATH"])
-
-geo_service  = GeoSvc("GeoSvc",
-        detectors=["{}/{}.xml".format(detector_path, detector_name)])
-podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=DEBUG)
-
-from Configurables import PodioInput
-from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
-from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
-from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
-
-from Configurables import Jug__Digi__ExampleCaloDigi as ExampleCaloDigi
-from Configurables import Jug__Digi__UFSDTrackerDigi as UFSDTrackerDigi
-from Configurables import Jug__Digi__EMCalorimeterDigi as EMCalorimeterDigi
-
-from Configurables import Jug__Base__MC2DummyParticle as MC2DummyParticle
-
-from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
-
-from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
-from Configurables import Jug__Reco__Tracker2SourceLinker as Tracker2SourceLinker
-#from Configurables import Jug__Reco__TrackerSourcesLinker as TrackerSourcesLinker
-#from Configurables import Jug__Reco__TrackingHitsSourceLinker as TrackingHitsSourceLinker
-from Configurables import Jug__Reco__TrackParamTruthInit as TrackParamTruthInit
-from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
-from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
-
-from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
-from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
-from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
-
-from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
-
-
-
-podioinput = PodioInput("PodioReader",
-        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","EcalBarrelHits"])#, OutputLevel=DEBUG)
-#"SiVertexBarrelHits",
-
-dummy = MC2DummyParticle("MC2Dummy",
-        inputCollection="mcparticles",
-        outputCollection="DummyReconstructedParticles",
-        smearing = 0.1)
-
-## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
-copier = MCCopier("MCCopier", 
-        inputCollection="mcparticles", 
-        outputCollection="mcparticles2") 
-trkcopier = TrkCopier("TrkCopier", 
-        inputCollection="TrackerBarrelHits", 
-        outputCollection="TrackerBarrelHits2") 
-
-ecal_digi = EMCalorimeterDigi("ecal_digi", 
-        inputHitCollection="EcalBarrelHits", 
-        outputHitCollection="RawEcalBarrelHits")
-
-ufsd_digi = UFSDTrackerDigi("ufsd_digi", 
-        inputHitCollection="TrackerBarrelHits",
-        outputHitCollection="TrackerBarrelRawHits",
-        timeResolution=8)
-ufsd_digi2 = UFSDTrackerDigi("ufsd_digi2", 
-        inputHitCollection="TrackerEndcapHits",
-        outputHitCollection="TrackerEndcapRawHits",
-        timeResolution=8)
-
-#vtx_digi = UFSDTrackerDigi("vtx_digi", 
-#        inputHitCollection="SiVertexBarrelHits",
-#        outputHitCollection="SiVertexBarrelRawHits",
-#        timeResolution=8)
-
-
-ecal_reco = EMCalReconstruction("ecal_reco", 
-        inputHitCollection=ecal_digi.outputHitCollection, 
-        outputHitCollection="RecEcalBarrelHits",
-        minModuleEdep=0.0*units.MeV,
-        OutputLevel=DEBUG)
-
-simple_cluster = SimpleClustering("simple_cluster", 
-        inputHitCollection=ecal_reco.outputHitCollection, 
-        outputClusterCollection="SimpleClusters",
-        outputProtoClusterCollection="SimpleProtoClusters",
-        minModuleEdep=1.0*units.MeV,
-        maxDistance=50.0*units.cm,
-        OutputLevel=DEBUG)
-
-trk_barrel_reco = TrackerHitReconstruction("trk_barrel_reco",
-        inputHitCollection="TrackerBarrelRawHits",
-        outputHitCollection="TrackerBarrelRecHits")
-
-trk_endcap_reco = TrackerHitReconstruction("trk_endcap_reco",
-        inputHitCollection="TrackerEndcapRawHits",
-        outputHitCollection="TrackerEndcapRecHits")
-
-#vtx_barrel_reco = TrackerHitReconstruction("vtx_barrel_reco",
-#        inputHitCollection = vtx_digi.outputHitCollection,
-#        outputHitCollection="VertexBarrelRecHits")
-
-# Source linker 
-sourcelinker = TrackerSourceLinker("trk_srclinker",
-        inputHitCollection="TrackerBarrelRecHits",
-        outputSourceLinks="BarrelTrackSourceLinks",
-        OutputLevel=DEBUG)
-
-trk_hits_srclnkr = Tracker2SourceLinker("trk_hits_srclnkr",
-        TrackerBarrelHits="TrackerBarrelRecHits",
-        TrackerEndcapHits="TrackerEndcapRecHits",
-        outputMeasurements="lnker2Measurements",
-        outputSourceLinks="lnker2Links",
-        allTrackerHits="linker2AllHits",
-        OutputLevel=DEBUG)
-
-## Track param init
-truth_trk_init = TrackParamTruthInit("truth_trk_init",
-        inputMCParticles="mcparticles",
-        outputInitialTrackParameters="InitTrackParams",
-        OutputLevel=DEBUG)
-
-clust_trk_init = TrackParamClusterInit("clust_trk_init",
-        inputClusters="SimpleClusters",
-        outputInitialTrackParameters="InitTrackParamsFromClusters",
-        OutputLevel=DEBUG)
-
-#vtxcluster_trk_init = TrackParamVertexClusterInit("vtxcluster_trk_init",
-#        inputVertexHits="VertexBarrelRecHits",
-#        inputClusters="SimpleClusters",
-#        outputInitialTrackParameters="InitTrackParamsFromVtxClusters",
-#        maxHitRadius=40.0*units.mm,
-#        OutputLevel=DEBUG)
-
-# Tracking algorithms
-trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
-        inputSourceLinks = sourcelinker.outputSourceLinks,
-        inputMeasurements = sourcelinker.outputMeasurements,
-        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
-        outputTrajectories="trajectories",
-        OutputLevel=DEBUG)
-parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
-        inputTrajectories="trajectories",
-        outputParticles="ReconstructedParticles",
-        outputTrackParameters="outputTrackParameters",
-        OutputLevel=DEBUG)
-
-trk_find_alg1 = TrackFindingAlgorithm("trk_find_alg1",
-        inputSourceLinks = trk_hits_srclnkr.outputSourceLinks,
-        inputMeasurements = trk_hits_srclnkr.outputMeasurements,
-        inputInitialTrackParameters= "InitTrackParamsFromClusters", 
-        outputTrajectories="trajectories1",
-        OutputLevel=DEBUG)
-parts_from_fit1 = ParticlesFromTrackFit("parts_from_fit1",
-        inputTrajectories="trajectories1",
-        outputParticles="ReconstructedParticles1",
-        outputTrackParameters="outputTrackParameters1",
-        OutputLevel=DEBUG)
-
-trk_find_alg2 = TrackFindingAlgorithm("trk_find_alg2",
-        inputSourceLinks = trk_hits_srclnkr.outputSourceLinks,
-        inputMeasurements = trk_hits_srclnkr.outputMeasurements,
-        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
-        #inputInitialTrackParameters= "InitTrackParamsFromVtxClusters", 
-        outputTrajectories="trajectories2",
-        OutputLevel=DEBUG)
-parts_from_fit2 = ParticlesFromTrackFit("parts_from_fit2",
-        inputTrajectories="trajectories2",
-        outputParticles="ReconstructedParticles2",
-        outputTrackParameters="outputTrackParameters2",
-        OutputLevel=DEBUG)
-
-
-#types = []
-## this printout is useful to check that the type information is passed to python correctly
-#print("---------------------------------------\n")
-#print("---\n# List of input and output types by class")
-#for configurable in sorted([ PodioInput, EICDataSvc, PodioOutput,
-#                             TrackerHitReconstruction,ExampleCaloDigi, 
-#                             UFSDTrackerDigi, TrackerSourceLinker,
-#                             PodioOutput],
-#                           key=lambda c: c.getType()):
-#    print("\"{}\":".format(configurable.getType()))
-#    props = configurable.getDefaultProperties()
-#    for propname, prop in sorted(props.items()):
-#        print(" prop name: {}".format(propname))
-#        if isinstance(prop, DataHandleBase):
-#            types.append(prop.type())
-#            print("  {}: \"{}\"".format(propname, prop.type()))
-#print("---")
-
-out = PodioOutput("out", filename=output_rec_file)
-out.outputCommands = ["keep *", 
-        "drop BarrelTrackSourceLinks", 
-        "drop InitTrackParams",
-        "drop trajectories",
-        "drop outputSourceLinks",
-        "drop outputInitialTrackParameters",
-        "drop mcparticles"
-        ]
-
-ApplicationMgr(
-    TopAlg = [podioinput, 
-              dummy,
-              copier, trkcopier,
-              ecal_digi, ufsd_digi2,ufsd_digi, #vtx_digi, 
-              ecal_reco, 
-              simple_cluster,
-              trk_barrel_reco, 
-              trk_endcap_reco, 
-              #vtx_barrel_reco, 
-              sourcelinker, trk_hits_srclnkr,
-              clust_trk_init, 
-              truth_trk_init, 
-              #vtxcluster_trk_init, 
-              trk_find_alg, parts_from_fit,
-              trk_find_alg1, parts_from_fit1,
-              trk_find_alg2, parts_from_fit2,
-              out
-              ],
-    EvtSel = 'NONE',
-    EvtMax   = n_events,
-    ExtSvc = [podioevent,geo_service],
-    OutputLevel=DEBUG
- )
-
-