diff --git a/CMakeLists.txt b/CMakeLists.txt
index bd8823e..3424b7c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,7 +2,7 @@
 # cmake file for building MarlinTrkProcessors
 # - based on $MARLIN/example/mymarlin/CMakeLists.txt
 #   by @author Jan Engels, Desy IT
-CMAKE_MINIMUM_REQUIRED(VERSION 2.6 FATAL_ERROR)
+CMAKE_MINIMUM_REQUIRED(VERSION 3.15 FATAL_ERROR)
 ########################################################
 
 
@@ -20,23 +20,6 @@ cmake_policy(SET CMP0008 NEW)
 
 FIND_PACKAGE( ILCUTIL REQUIRED COMPONENTS ILCSOFT_CMAKE_MODULES )
 
-# load default settings from ILCSOFT_CMAKE_MODULES
-INCLUDE( ilcsoft_default_settings )
-
-SET( COMPILER_FLAGS "-Wno-effc++" )
-MESSAGE( STATUS "FLAGS ${COMPILER_FLAGS}" )
-FOREACH( FLAG ${COMPILER_FLAGS} )
-  CHECK_CXX_COMPILER_FLAG( "${FLAG}" CXX_FLAG_WORKS_${FLAG} )
-  IF( ${CXX_FLAG_WORKS_${FLAG}} )
-    MESSAGE ( STATUS "Adding ${FLAG} to CXX_FLAGS" )
-    ### We append the flag, so they overwrite things from somewhere else
-    SET ( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAG} ")
-  ELSE()
-    MESSAGE ( STATUS "NOT Adding ${FLAG} to CXX_FLAGS" )
-  ENDIF()
-ENDFOREACH()
-
-
 FIND_PACKAGE( Marlin 1.0 REQUIRED ) # minimum required Marlin version
 INCLUDE_DIRECTORIES( SYSTEM ${Marlin_INCLUDE_DIRS} )
 LINK_LIBRARIES( ${Marlin_LIBRARIES} )
@@ -128,16 +111,5 @@ INCLUDE_DIRECTORIES( . ./source/Digitisers/include ./source/Refitting/include ./
 # add library
 FILE( GLOB_RECURSE library_sources "source/*.cc" )
 
-# needed for adding header files to xcode project
-IF(CMAKE_GENERATOR MATCHES "Xcode")
-  FILE( GLOB_RECURSE library_headers "*.h" )
-  ADD_SHARED_LIBRARY( ${PROJECT_NAME} ${library_sources} ${library_headers})
-ELSE()
-  ADD_SHARED_LIBRARY( ${PROJECT_NAME} ${library_sources} )
-ENDIF()  
-
-INSTALL_SHARED_LIBRARY( ${PROJECT_NAME} DESTINATION lib )
-
-# display some variables and write them to cache
-DISPLAY_STD_VARIABLES()
+ADD_LIBRARY( ${PROJECT_NAME} SHARED ${library_sources} )
 
diff --git a/source/Utils/include/FilterClusters.h b/source/Utils/include/FilterClusters.h
index ab76767..7d6af94 100644
--- a/source/Utils/include/FilterClusters.h
+++ b/source/Utils/include/FilterClusters.h
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <marlin/Processor.h>
+#include <TH1F.h>
 
 //#include <ACTSTracking/GeometryIdMappingTool.hxx>
 
@@ -39,18 +40,33 @@ class FilterClusters : public marlin::Processor
    //! Input track collection
    std::string _InTrackerHitCollection {};
    std::string _InRelationCollection {};
+   std::string _InSimTrackerHitCollection {};
 
    //! Output track collection
    std::string _OutTrackerHitCollection {};
    std::string _OutRelationCollection {};
+   std::string _OutSimTrackerHitCollection {};
 
    //! Ranges for theta
-   std::vector<std::string> _ThetaRanges;
+  std::vector<std::string> _ThetaRanges;
 
    //! Cut-offs for cluster size in various theta ranges
-   std::vector<std::string> _ClusterSize;
-
-   //! Layers to be filtered
-   std::vector<std::string> _Layers;
+  std::vector<std::string> _ClusterSize;
+
+  //! Layers to be filtered
+  std::vector<std::string> _Layers;
+  
+  //!Number of bins in theta
+  std::string _ThetaBins;
+
+  bool m_fillHistos{};
+
+  // --- Diagnostic histograms:
+  TH1F *m_clusterTheta_beforeCut = nullptr;
+  TH1F *m_clusterTheta_afterCut = nullptr;
+  TH1F *m_clusterLayer_beforeCut = nullptr;
+  TH1F *m_clusterLayer_afterCut = nullptr;
+  TH1F *m_clusterSize_beforeCut = nullptr;
+  TH1F *m_clusterSize_afterCut = nullptr;
 
 };
diff --git a/source/Utils/include/FilterTimeHits.h b/source/Utils/include/FilterTimeHits.h
index d98140a..871d8d6 100644
--- a/source/Utils/include/FilterTimeHits.h
+++ b/source/Utils/include/FilterTimeHits.h
@@ -6,6 +6,9 @@
 #include <string>
 #include <vector>
 
+#include <EVENT/TrackerHitPlane.h>
+#include <EVENT/SimTrackerHit.h>
+
 #include "DDRec/Surface.h"
 #include "DDRec/SurfaceManager.h"
 
@@ -55,10 +58,12 @@ class FilterTimeHits : public Processor
 
 protected:
     // --- Input/output collection names:
-    std::vector<std::string> m_inputTrackerHitsCollNames{};
-    std::vector<std::string> m_inputTrackerSimHitsCollNames{};
-    std::vector<std::string> m_inputTrackerHitRelNames{};
+    std::vector<std::string> m_inputTrackerHitsCollNames{}; /// Reconstructed clusters
+    std::vector<std::string> m_inputTrackerHitsConstituentsCollNames{}; // Individual hits in clusters (if available)
+    std::vector<std::string> m_inputTrackerSimHitsCollNames{}; // Simulated energy deposits (if available)
+    std::vector<std::string> m_inputTrackerHitRelNames{}; // Relations between reconstructed clusters and simulated ones (if available)
     std::vector<std::string> m_outputTrackerHitsCollNames{};
+    std::vector<std::string> m_outputTrackerHitsConstituentsCollNames{};
     std::vector<std::string> m_outputTrackerSimHitsCollNames{};
     std::vector<std::string> m_outputTrackerHitRelNames{};
 
@@ -69,7 +74,14 @@ class FilterTimeHits : public Processor
     double m_time_max{};
 
     // --- Diagnostic histograms:
-    TH1F *m_corrected_time = nullptr;
+    TH1F *m_corrected_time_before = nullptr;
+    TH1F *m_corrected_time_after = nullptr;
+
+    // ---- Utility functions
+    // Make a full copy of a TrackerHitPlane object (no copy-constructor defined)
+    TrackerHitPlane* copyTrackerHitPlane(TrackerHitPlane *hit);
+    // Make a fully copy of a SimTrackerHit object
+    SimTrackerHit* copySimTrackerHit(SimTrackerHit *hit);
 
     // --- Run and event counters:
     int _nRun{};
diff --git a/source/Utils/src/FilterClusters.cc b/source/Utils/src/FilterClusters.cc
index 1797855..c3ee2ab 100644
--- a/source/Utils/src/FilterClusters.cc
+++ b/source/Utils/src/FilterClusters.cc
@@ -6,11 +6,14 @@
 
 #include <EVENT/Track.h>
 #include <EVENT/TrackerHit.h>
+#include <IMPL/TrackerHitPlaneImpl.h>
 #include <EVENT/SimTrackerHit.h>
+#include <IMPL/SimTrackerHitImpl.h>
 #include <IMPL/LCCollectionVec.h>
 #include <UTIL/LCTrackerConf.h>
 #include <UTIL/CellIDDecoder.h>
 #include <IMPL/LCRelationImpl.h>
+#include <IMPL/LCFlagImpl.h>
 
 #include <marlin/AIDAProcessor.h>
 
@@ -18,167 +21,347 @@
 
 #include "marlin/VerbosityLevels.h"
 
-
-FilterClusters aFilterClusters ;
+FilterClusters aFilterClusters;
 
 FilterClusters::FilterClusters()
-  : Processor("FilterClusters")
+    : Processor("FilterClusters")
 {
   // modify processor description
   _description = "FilterClusters processor filters a collection of tracker hits based on cluster size and outputs a filtered collection";
 
   // register steering parameters: name, description, class-variable, default value
   registerProcessorParameter("ThetaRanges",
-         "Divide theta into bins for different cluster size cuts",
-         _ThetaRanges, 
-         {}
-          );
+                             "Divide theta into bins for different cluster size cuts",
+                             _ThetaRanges,
+                             {});
+  registerProcessorParameter("ThetaBins",
+                             "Number of bins in theta",
+                             _ThetaBins,
+                             {});
   registerProcessorParameter("ClusterSize",
-		  	 "Maximum cluster size for each theta range",
-			   _ClusterSize,
-			   {}
-			    );  
+                             "Maximum cluster size for each theta range",
+                             _ClusterSize,
+                             {});
   registerProcessorParameter("Layers",
-		  	 "Layers to be filtered",
-			   _Layers,
-			   {}
-			    );  
-  registerInputCollection( LCIO::TRACKERHIT,
-		  	 "InTrackerHitCollection" ,
-			   "Name of the input tracker hit collection",
-			   _InTrackerHitCollection,
-		     _InTrackerHitCollection
-		 	    );
-
-  registerInputCollection( LCIO::LCRELATION,
-		     "InRelationCollection" ,
-			   "Name of the input relation collection",
-			   _InRelationCollection,
-		     _InRelationCollection
-		 	    );
-
-  registerOutputCollection( LCIO::TRACKERHIT,
-		  	 "OutTrackerHitCollection" ,
-			   "Name of output tracker hit collection",
-			   _OutTrackerHitCollection,
-			   std::string("FilteredVBTrackerHits")
-			    );
-
-    registerOutputCollection( LCIO::LCRELATION,
-		  	 "OutRelationCollection" ,
-			   "Name of output relation collection",
-			   _OutRelationCollection,
-			   std::string("FilteredVBTrackerHitsRelations")
-			    );
+                             "Layers to be filtered",
+                             _Layers,
+                             {});
+  registerInputCollection(LCIO::SIMTRACKERHIT,
+                          "InSimTrackerHitCollection",
+                          "Name of the input sim tracker hit collection",
+                          _InSimTrackerHitCollection,
+                          _InSimTrackerHitCollection);
+
+  registerInputCollection(LCIO::TRACKERHITPLANE,
+                          "InTrackerHitCollection",
+                          "Name of the input tracker hit collection",
+                          _InTrackerHitCollection,
+                          _InTrackerHitCollection);
+
+  registerInputCollection(LCIO::LCRELATION,
+                          "InRelationCollection",
+                          "Name of the input relation collection",
+                          _InRelationCollection,
+                          _InRelationCollection);
+
+  registerOutputCollection(LCIO::SIMTRACKERHIT,
+                           "OutSimTrackerHitCollection",
+                           "Name of output sim tracker hit collection",
+                           _OutSimTrackerHitCollection,
+                           _OutSimTrackerHitCollection);
+
+  registerOutputCollection(LCIO::TRACKERHITPLANE,
+                           "OutTrackerHitCollection",
+                           "Name of output tracker hit collection",
+                           _OutTrackerHitCollection,
+                           _OutTrackerHitCollection);
 
+  registerOutputCollection(LCIO::LCRELATION,
+                           "OutRelationCollection",
+                           "Name of output relation collection",
+                           _OutRelationCollection,
+                           _OutRelationCollection);
+
+  registerProcessorParameter("FillHistograms",
+                             "Flag to fill the diagnostic histograms",
+                             m_fillHistos,
+                             false);
 }
 
 void FilterClusters::init()
 {
   // Print the initial parameters
-  printParameters() ;
-}
+  printParameters();
 
-void FilterClusters::processRunHeader( LCRunHeader* /*run*/)
-{ }
+  marlin::AIDAProcessor::histogramFactory(this);
+
+  m_clusterTheta_beforeCut = new TH1F("m_ClusterTheta_before", "Cluster Theta [radian]", 32, 0., 3.2);
+  m_clusterTheta_afterCut = new TH1F("m_ClusterTheta_after", "Cluster Theta [radian]", 32, 0., 3.2);
+  m_clusterLayer_beforeCut = new TH1F("m_ClusterLayer_before", "Cluster Layer Index", 10, 0, 10);
+  m_clusterLayer_afterCut = new TH1F("m_ClusterLayer_after", "Cluster Layer Index", 10, 0, 10);
+  m_clusterSize_beforeCut = new TH1F("m_ClusterSize_before", "Cluster hit multiplicity", 20, 0, 20);
+  m_clusterSize_afterCut = new TH1F("m_ClusterSize_after", "Cluster hit multiplicity", 20, 0, 20);
+}
 
-void FilterClusters::processEvent( LCEvent * evt )
+void FilterClusters::processRunHeader(LCRunHeader * /*run*/)
 {
-  // Make the output track collection
-  LCCollectionVec *OutTrackerHitCollection = new LCCollectionVec(LCIO::TRACKERHIT);
-  OutTrackerHitCollection->setSubset(true);
-  LCCollectionVec *OutRelationCollection   = new LCCollectionVec(LCIO::LCRELATION);
-  OutRelationCollection->setSubset(true);
+}
 
+void FilterClusters::processEvent(LCEvent *evt)
+{
   // Get input collection
-  LCCollection* InTrackerHitCollection  = evt->getCollection(_InTrackerHitCollection);
-  LCCollection* InRelationCollection    = evt->getCollection(_InRelationCollection);
+  LCCollection *InTrackerHitCollection = evt->getCollection(_InTrackerHitCollection);
+  LCCollection *InRelationCollection = evt->getCollection(_InRelationCollection);
+  LCCollection *InSimTrackerHitCollection = evt->getCollection(_InSimTrackerHitCollection);
 
-  if( InTrackerHitCollection->getTypeName() != lcio::LCIO::TRACKERHITPLANE )
-    { throw EVENT::Exception( "Invalid collection type: " + InTrackerHitCollection->getTypeName() ) ; }
+  streamlog_out(DEBUG0) << "Starting processing event." << std::endl;
+
+  if (InTrackerHitCollection->getTypeName() != lcio::LCIO::TRACKERHITPLANE)
+  {
+    throw EVENT::Exception("Invalid collection type: " + InTrackerHitCollection->getTypeName());
     streamlog_out(DEBUG0) << "Wrong collection type for TrackerHitCollection. \n";
-  if( InRelationCollection->getTypeName() != lcio::LCIO::LCRELATION )
-    { throw EVENT::Exception( "Invalid collection type: " + InRelationCollection->getTypeName() ) ; }
+  }
+  if (InRelationCollection->getTypeName() != lcio::LCIO::LCRELATION)
+  {
+    throw EVENT::Exception("Invalid collection type: " + InRelationCollection->getTypeName());
     streamlog_out(DEBUG0) << "Wrong collection type for InRelationCollection. \n";
+  }
+  if (InSimTrackerHitCollection->getTypeName() != lcio::LCIO::SIMTRACKERHIT)
+  {
+    throw EVENT::Exception("Invalid collection type: " + InSimTrackerHitCollection->getTypeName());
+    streamlog_out(DEBUG0) << "Wrong collection type for SimTrackerHitCollection. \n";
+  }
+
+  streamlog_out(DEBUG1) << "Number of Elements in Tracker Hits Collection: " << InTrackerHitCollection->getNumberOfElements() << std::endl;
 
+  // Make the output collections: reco hits, sim hits, reco-sim relationship
+  std::string encoderString = InTrackerHitCollection->getParameters().getStringVal("CellIDEncoding");
+  LCCollectionVec *OutTrackerHitCollection = new LCCollectionVec(LCIO::TRACKERHITPLANE);
+  OutTrackerHitCollection->parameters().setValue("CellIDEncoding", encoderString);
+  LCFlagImpl lcFlag(InTrackerHitCollection->getFlag());
+  OutTrackerHitCollection->setFlag(lcFlag.getFlag());
+  // OutTrackerHitCollection->setSubset(true);
+
+  LCCollectionVec *OutSimTrackerHitCollection = new LCCollectionVec(LCIO::SIMTRACKERHIT);
+  OutSimTrackerHitCollection->parameters().setValue("CellIDEncoding", encoderString);
+  LCFlagImpl lcFlag_sim(InSimTrackerHitCollection->getFlag());
+  OutSimTrackerHitCollection->setFlag(lcFlag_sim.getFlag());
+  // OutSimTrackerHitCollection->setSubset(true);
+
+  LCCollectionVec *OutRelationCollection = new LCCollectionVec(LCIO::LCRELATION);
+  LCFlagImpl lcFlag_rel(InRelationCollection->getFlag());
+  OutRelationCollection->setFlag(lcFlag_rel.getFlag());
+  // OutRelationCollection->setSubset(true);
 
-  streamlog_out(DEBUG0) << "Number of Elements in VB Tracker Hits Collection: " << InTrackerHitCollection->getNumberOfElements() <<std::endl;
   // Filter
-  for(int i=0; i<InTrackerHitCollection->getNumberOfElements(); ++i) //loop through all hits
+  for (int i = 0; i < InTrackerHitCollection->getNumberOfElements(); ++i) // loop through all hits
+  {
+    streamlog_out(DEBUG2) << "Loop over hits:" << i << std::endl;
+    TrackerHitPlane *trkhit = static_cast<TrackerHitPlane *>(InTrackerHitCollection->getElementAt(i)); // define trkhit var, pointer to i'th element of tracker hits
+
+    if (!trkhit)
     {
-      streamlog_out(DEBUG0) << "Loop over hits opened. \n";
-      EVENT::TrackerHit *trkhit=static_cast<EVENT::TrackerHit*>(InTrackerHitCollection->getElementAt(i)); //define trkhit var, pointer to i'th element of tracker hits
-      EVENT::LCRelation *rel=static_cast<EVENT::LCRelation*>(InRelationCollection->getElementAt(i));
-      
-      //Calculating theta 
-      float x = trkhit->getPosition()[0];
-      float y = trkhit->getPosition()[1];
-      float z = trkhit->getPosition()[2];
-      float r = sqrt(pow(x,2)+pow(y,2));
-      float incidentTheta = std::atan(r/z); 
-
-      if(incidentTheta<0)
-        incidentTheta += M_PI;
-
-      //Calculating cluster size
-      const lcio::LCObjectVec &rawHits = trkhit->getRawHits();
-      float max = -1000000;
-      float min = 1000000;
-      for (size_t j=0; j<rawHits.size(); ++j) {
-        lcio::SimTrackerHit *hitConstituent = dynamic_cast<lcio::SimTrackerHit*>( rawHits[j] );
-        const double *localPos = hitConstituent->getPosition();
-        x = localPos[0];
-        y = localPos[1];
-        if (y < min){
-          min = y;
-        }
-        else if (y > max){
-          max = y;          
-        }
+      streamlog_out(WARNING) << "Cannot retrieve valid point to cluster. Skipping it" << std::endl;
+    }
+
+    // Calculating theta
+    float x = trkhit->getPosition()[0];
+    float y = trkhit->getPosition()[1];
+    float z = trkhit->getPosition()[2];
+    float r = sqrt(pow(x, 2) + pow(y, 2));
+    float incidentTheta = std::atan(r / z);
+
+    if (incidentTheta < 0)
+      incidentTheta += M_PI;
+
+    // Calculating cluster size
+    const lcio::LCObjectVec &rawHits = trkhit->getRawHits();
+    float ymax = -1000000;
+    float xmax = -1000000;
+    float ymin = 1000000;
+    float xmin = 1000000;
+    streamlog_out(DEBUG1) << "Looping over hits constituents." << std::endl;
+    for (size_t j = 0; j < rawHits.size(); ++j)
+    {
+      lcio::SimTrackerHit *hitConstituent = dynamic_cast<lcio::SimTrackerHit *>(rawHits[j]);
+      const double *localPos = hitConstituent->getPosition();
+      float x_local = localPos[0];
+      float y_local = localPos[1];
+
+      if (y_local < ymin)
+      {
+        ymin = y_local;
       }
-      float cluster_size = (max - min)+1;
-
-      //Get hit subdetector/layer 
-      std::string _encoderString = lcio::LCTrackerCellID::encoding_string();
-      UTIL::CellIDDecoder<lcio::TrackerHit> decoder(_encoderString);
-      uint32_t systemID = decoder(trkhit)["system"];
-      uint32_t layerID = decoder(trkhit)["layer"];
-      bool filter_layer = false;
-      for (size_t j=0; j<_Layers.size(); ++j){
-        if (layerID == std::stof(_Layers[j])) {
-          filter_layer = true;
-        }
+      if (y_local > ymax)
+      {
+        ymax = y_local;
+      }
+
+      if (x_local < xmin)
+      {
+        xmin = x_local;
+      }
+      if (x_local > xmax)
+      {
+        xmax = x_local;
+      }
+    }
+    float cluster_size_y = (ymax - ymin) + 1;
+    float cluster_size_x = (xmax - xmin) + 1;
+    float cluster_size = rawHits.size();
+
+    streamlog_out(DEBUG2) << "Cluster size:" << cluster_size << std::endl;
+
+    // Get hit subdetector/layer
+    std::string _encoderString = lcio::LCTrackerCellID::encoding_string();
+    UTIL::CellIDDecoder<lcio::TrackerHit> decoder(_encoderString);
+    uint32_t systemID = decoder(trkhit)["system"];
+    uint32_t layerID = decoder(trkhit)["layer"];
+    bool filter_layer = false;
+
+    streamlog_out(DEBUG3) << "Decoded system/layer:" << systemID << "/" << layerID << std::endl;
+
+    int rows = _Layers.size(), cols = std::stoi(_ThetaBins);
+    if ((rows * (cols + 1) != _ThetaRanges.size()) || (rows * cols != _ClusterSize.size()))
+    {
+      std::cout << "Either theta cuts or cluster cuts not provided for each layer. Please change the config, exiting now..." << std::endl;
+      return;
+    }
+
+    std::vector<std::vector<float>> _thetaBins_byLayer;
+    std::vector<std::vector<float>> _clusterSizeCuts_byLayer;
+
+    for (int k = 0; k < rows; ++k)
+    {
+      std::vector<float> row;
+      for (int j = 0; j <= cols; ++j)
+      {
+        row.push_back(std::stof(_ThetaRanges[j + k*(cols+1)]));
+      }
+      _thetaBins_byLayer.push_back(row);
+    }
+
+    for (int k = 0; k < rows; ++k)
+    {
+      std::vector<float> row;
+      for (int j = 0; j < cols; ++j)
+      {
+        row.push_back(std::stof(_ClusterSize[j + k*cols]));
       }
-      streamlog_out(DEBUG0) << "Filter layer: " << filter_layer << std::endl;
-      
-      for (size_t j=0; j<_ThetaRanges.size()-1; ++j) {
-        streamlog_out( DEBUG0 ) << "theta: " << incidentTheta << std::endl;
-        float min_theta = std::stof(_ThetaRanges[j]);
-        float max_theta = std::stof(_ThetaRanges[j+1]);
-        streamlog_out( DEBUG0 ) << "theta range: " << min_theta << ", " << max_theta << std::endl;
-
-        if(incidentTheta > min and incidentTheta <= max and not filter_layer){
-          streamlog_out( DEBUG0 ) << "theta in range" << std::endl;
-          streamlog_out( DEBUG0 ) << "cluster size cut off: " << _ClusterSize[j] << std::endl;
-          streamlog_out( DEBUG0 ) << "cluster size: " << cluster_size << std::endl;
-          if(cluster_size < std::stof(_ClusterSize[j])) {
-            streamlog_out( DEBUG0 ) << "cluster added" << std::endl;
-            OutTrackerHitCollection->addElement(trkhit); 
-            OutRelationCollection->addElement(rel); }
-          else {
-            streamlog_out( DEBUG0 ) << "cluster rejected" << std::endl;
+      _clusterSizeCuts_byLayer.push_back(row);
+    }
+
+    int layerInd = -1;
+    for (size_t j = 0; j < _Layers.size(); ++j)
+    {
+      if (layerID == std::stof(_Layers[j]))
+      {
+        filter_layer = true;
+	layerInd = j;
+        break;
+      }
+    }
+    streamlog_out(DEBUG1) << "Filter layer: " << filter_layer << std::endl;
+
+    if (m_fillHistos)
+    {
+      m_clusterTheta_beforeCut->Fill(incidentTheta);
+      m_clusterLayer_beforeCut->Fill(layerID);
+      m_clusterSize_beforeCut->Fill(cluster_size);
+    }
+
+    bool store_hit = true;
+    if (filter_layer)
+    {
+      store_hit = false;
+      for (size_t j = 0; j < _thetaBins_byLayer[layerInd].size() - 1; ++j)
+      {
+        streamlog_out(DEBUG0) << "theta: " << incidentTheta << std::endl;
+        float min_theta = _thetaBins_byLayer[layerInd][j];
+        float max_theta = _thetaBins_byLayer[layerInd][j + 1];
+        streamlog_out(DEBUG0) << "theta range: " << min_theta << ", " << max_theta << std::endl;
+
+        if (incidentTheta >= min_theta and incidentTheta <= max_theta and filter_layer)
+        {
+          streamlog_out(DEBUG0) << "theta in range" << std::endl;
+          streamlog_out(DEBUG0) << "cluster size cut off: " << _clusterSizeCuts_byLayer[layerInd][j] << std::endl;
+          streamlog_out(DEBUG0) << "cluster size: " << cluster_size << std::endl;
+          if (cluster_size < _clusterSizeCuts_byLayer[layerInd][j])
+          {
+	    store_hit = true;
+            streamlog_out(DEBUG0) << "Adding reco/sim clusters and relation to output collections" << std::endl;
+	    break;
           }
         }
-        else{
-          streamlog_out( DEBUG0 ) << "theta out of range or layer filtered" << std::endl;
-        }
       }
     }
 
+    if (store_hit)
+    {
+      if (m_fillHistos)
+      {
+        m_clusterTheta_afterCut->Fill(incidentTheta);
+        m_clusterLayer_afterCut->Fill(layerID);
+        m_clusterSize_afterCut->Fill(cluster_size);
+      }
+
+      EVENT::LCRelation *rel = static_cast<EVENT::LCRelation *>(InRelationCollection->getElementAt(i));
+      SimTrackerHit *simhit = dynamic_cast<SimTrackerHit *>(rel->getTo());
+      SimTrackerHitImpl *simhit_new = new SimTrackerHitImpl();
+      simhit_new->setCellID0(simhit->getCellID0());
+      simhit_new->setCellID1(simhit->getCellID1());
+      simhit_new->setPosition(simhit->getPosition());
+      simhit_new->setEDep(simhit->getEDep());
+      simhit_new->setTime(simhit->getTime());
+      simhit_new->setMCParticle(simhit->getMCParticle());
+      simhit_new->setMomentum(simhit->getMomentum());
+      simhit_new->setPathLength(simhit->getPathLength());
+      simhit_new->setQuality(simhit->getQuality());
+      simhit_new->setOverlay(simhit->isOverlay());
+      simhit_new->setProducedBySecondary(simhit->isProducedBySecondary());
+      OutSimTrackerHitCollection->addElement(simhit_new);
+
+      TrackerHitPlaneImpl *hit_new = new TrackerHitPlaneImpl();
+      hit_new->setCellID0(trkhit->getCellID0());
+      hit_new->setCellID1(trkhit->getCellID1());
+      hit_new->setType(trkhit->getType());
+      hit_new->setPosition(trkhit->getPosition());
+      hit_new->setU(trkhit->getU());
+      hit_new->setV(trkhit->getV());
+      hit_new->setdU(trkhit->getdU());
+      hit_new->setdV(trkhit->getdV());
+      hit_new->setEDep(trkhit->getEDep());
+      hit_new->setEDepError(trkhit->getEDepError());
+      hit_new->setTime(trkhit->getTime());
+      hit_new->setQuality(trkhit->getQuality());
+      const lcio::LCObjectVec &rawHits_new = trkhit->getRawHits();
+      for (size_t k = 0; k < rawHits_new.size(); ++k)
+      {
+        lcio::SimTrackerHit *hitConstituent = dynamic_cast<lcio::SimTrackerHit *>(rawHits_new[k]);
+        hit_new->rawHits().push_back(hitConstituent);
+      }
+      OutTrackerHitCollection->addElement(hit_new);
+
+      LCRelationImpl *rel_new = new LCRelationImpl();
+      rel_new->setFrom(hit_new);
+      rel_new->setTo(simhit_new);
+      rel_new->setWeight(rel->getWeight());
+      OutRelationCollection->addElement(rel_new);
+    }
+    else
+    {
+      streamlog_out(DEBUG0) << "cluster rejected" << std::endl;
+    }
+  }
+
   // Save output track collection
-  evt->addCollection(OutTrackerHitCollection, _OutTrackerHitCollection); 
-  evt->addCollection(OutRelationCollection, _OutRelationCollection); 
+  evt->addCollection(OutTrackerHitCollection, _OutTrackerHitCollection);
+  evt->addCollection(OutRelationCollection, _OutRelationCollection);
+  evt->addCollection(OutSimTrackerHitCollection, _OutSimTrackerHitCollection);
+
+  streamlog_out(MESSAGE) << " " << OutTrackerHitCollection->size() << " reco clusters added to the collection: " << _OutTrackerHitCollection << ", " << OutSimTrackerHitCollection->size() << " sim hits added to the collection: " << _OutSimTrackerHitCollection << " and " << OutRelationCollection->size() << " relations added to the collection: " << _OutRelationCollection << std::endl;
 }
 
 void FilterClusters::end()
-{ }
+{
+}
diff --git a/source/Utils/src/FilterTimeHits.cc b/source/Utils/src/FilterTimeHits.cc
index 1419668..2f95f37 100644
--- a/source/Utils/src/FilterTimeHits.cc
+++ b/source/Utils/src/FilterTimeHits.cc
@@ -32,6 +32,11 @@ FilterTimeHits::FilterTimeHits() : Processor("FilterTimeHits")
                                m_inputTrackerHitsCollNames,
                                {});
 
+    registerProcessorParameter("TrackerHitConstituentsInputCollections",
+                               "Name of the tracker hit constituents input collections",
+                               m_inputTrackerHitsConstituentsCollNames,
+                               {});
+
     registerProcessorParameter("TrackerSimHitInputCollections",
                                "Name of the tracker simhit input collections",
                                m_inputTrackerSimHitsCollNames,
@@ -47,6 +52,11 @@ FilterTimeHits::FilterTimeHits() : Processor("FilterTimeHits")
                                m_outputTrackerHitsCollNames,
                                {});
 
+    registerProcessorParameter("TrackerHitConstituentsOutputCollections",
+                               "Name of the tracker hit output collections",
+                               m_outputTrackerHitsConstituentsCollNames,
+                               {});
+
     registerProcessorParameter("TrackerSimHitOutputCollections",
                                "Name of the tracker simhit output collections",
                                m_outputTrackerSimHitsCollNames,
@@ -96,7 +106,8 @@ void FilterTimeHits::init()
 
     AIDAProcessor::histogramFactory(this);
 
-    m_corrected_time = new TH1F("m_corrected_time", "Corrected time of the hit [ps]", 1000, -250., 250.);
+    m_corrected_time_before = new TH1F("m_corrected_time_before", "Corrected time of the hit before filter [ns]", 1000, -25., 25.);
+    m_corrected_time_after = new TH1F("m_corrected_time_after", "Corrected time of the hit after filter [ns]", 1000, -25., 25.);
 }
 
 void FilterTimeHits::processRunHeader(LCRunHeader *run)
@@ -104,6 +115,72 @@ void FilterTimeHits::processRunHeader(LCRunHeader *run)
     _nRun++;
 }
 
+TrackerHitPlane *FilterTimeHits::copyTrackerHitPlane(TrackerHitPlane *hit)
+{
+    // Create new object
+    TrackerHitPlaneImpl *hit_new = new TrackerHitPlaneImpl();
+
+    // Make a copy (no copy-constructor)
+    hit_new->setCellID0(hit->getCellID0());
+    hit_new->setCellID1(hit->getCellID1());
+    hit_new->setType(hit->getType());
+    hit_new->setPosition(hit->getPosition());
+    hit_new->setU(hit->getU());
+    hit_new->setV(hit->getV());
+    hit_new->setdU(hit->getdU());
+    hit_new->setdV(hit->getdV());
+    hit_new->setEDep(hit->getEDep());
+    hit_new->setEDepError(hit->getEDepError());
+    hit_new->setTime(hit->getTime());
+    hit_new->setQuality(hit->getQuality());
+    const_cast<EVENT::FloatVec &>(hit_new->getCovMatrix()) = hit->getCovMatrix(); // no setter for covariance matrix?
+
+    // need to clone individual hits, if present
+    const lcio::LCObjectVec &rawHits = hit->getRawHits();
+    for (size_t j = 0; j < rawHits.size(); ++j)
+    {
+        // Use (default) copy-constructor of SimTrackerHitImpl
+        lcio::SimTrackerHit *hitConstituent = copySimTrackerHit(dynamic_cast<SimTrackerHit *>(rawHits[j]));
+        if (!hitConstituent)
+        {
+            static bool first_error = true;
+            if (first_error)
+            {
+                streamlog_out(WARNING) << "Cannot access individual hits of reco ID clusters. Skipping (no further WARNING will be printed.)" << std::endl;
+                continue;
+            }
+        }
+        hit_new->rawHits().push_back(hitConstituent);
+    }
+
+    // return new object
+    return hit_new;
+}
+
+SimTrackerHit *FilterTimeHits::copySimTrackerHit(SimTrackerHit *hit)
+{
+    // use (default) copy-constructor
+    lcio::SimTrackerHitImpl *simhit_new = new lcio::SimTrackerHitImpl(*(dynamic_cast<lcio::SimTrackerHitImpl *>(hit)));
+
+    /*
+    // manual copy, deprecated
+    SimTrackerHitImpl *simhit_new = new SimTrackerHitImpl();
+    simhit_new->setCellID0(simhit->getCellID0());
+    simhit_new->setCellID1(simhit->getCellID1());
+    simhit_new->setPosition(simhit->getPosition());
+    simhit_new->setEDep(simhit->getEDep());
+    simhit_new->setTime(simhit->getTime());
+    simhit_new->setMCParticle(simhit->getMCParticle());
+    simhit_new->setMomentum(simhit->getMomentum());
+    simhit_new->setPathLength(simhit->getPathLength());
+    simhit_new->setQuality(simhit->getQuality());
+    simhit_new->setOverlay(simhit->isOverlay());
+    simhit_new->setProducedBySecondary(simhit->isProducedBySecondary());
+    */
+
+    return simhit_new;
+}
+
 void FilterTimeHits::processEvent(LCEvent *evt)
 {
 
@@ -140,10 +217,12 @@ void FilterTimeHits::processEvent(LCEvent *evt)
 
     const unsigned int nTrackerHitCol = m_inputTrackerHitsCollNames.size();
     std::vector<LCCollection *> inputHitColls(nTrackerHitCol);
+    std::vector<LCCollection *> inputHitConstituentsColls(nTrackerHitCol);
     std::vector<LCCollection *> inputSimHitColls(nTrackerHitCol);
     std::vector<LCCollection *> inputHitRels(nTrackerHitCol);
 
     std::vector<LCCollectionVec *> outputTrackerHitColls(nTrackerHitCol);
+    std::vector<LCCollectionVec *> outputTrackerHitConstituentsColls(nTrackerHitCol);
     std::vector<LCCollectionVec *> outputTrackerSimHitColls(nTrackerHitCol);
     std::vector<LCCollectionVec *> outputTrackerHitRels(nTrackerHitCol);
 
@@ -159,31 +238,56 @@ void FilterTimeHits::processEvent(LCEvent *evt)
         {
             streamlog_out(WARNING) << m_inputTrackerHitsCollNames[icol]
                                    << " collection not available" << std::endl;
-            continue;
+            continue; // this is mandatory to do anything
+        }
+
+        // get the reco hits contituents (if available)
+        bool hasHitConstituents = false;
+        try
+        {
+            if (m_inputTrackerHitsConstituentsCollNames.size() > 0)
+                if (m_inputTrackerHitsConstituentsCollNames[icol] != "")
+                    hasHitConstituents = true;
+            if (hasHitConstituents)
+                inputHitConstituentsColls[icol] = evt->getCollection(m_inputTrackerHitsConstituentsCollNames[icol]);
+            else
+                inputHitConstituentsColls[icol] = nullptr;
+        }
+        catch (lcio::DataNotAvailableException &e)
+        {
+            streamlog_out(WARNING) << m_inputTrackerHitsConstituentsCollNames[icol]
+                                   << " collection not available. This is expected in simplified digitization settings." << std::endl;
+            inputHitConstituentsColls[icol] = nullptr; // optional collection
         }
 
         // get the sim hits
         try
         {
-            inputSimHitColls[icol] = evt->getCollection(m_inputTrackerSimHitsCollNames[icol]);
+            if (m_inputTrackerSimHitsCollNames[icol] != "")
+                inputSimHitColls[icol] = evt->getCollection(m_inputTrackerSimHitsCollNames[icol]);
+            else
+                inputSimHitColls[icol] = nullptr;
         }
         catch (lcio::DataNotAvailableException &e)
         {
             streamlog_out(WARNING) << m_inputTrackerSimHitsCollNames[icol]
                                    << " collection not available" << std::endl;
-            continue;
+            inputSimHitColls[icol] = nullptr; // optional collection
         }
 
         // get the reco-sim relations
         try
         {
-            inputHitRels[icol] = evt->getCollection(m_inputTrackerHitRelNames[icol]);
+            if (m_inputTrackerHitRelNames[icol] != "")
+                inputHitRels[icol] = evt->getCollection(m_inputTrackerHitRelNames[icol]);
+            else
+                inputHitRels[icol] = nullptr;
         }
         catch (lcio::DataNotAvailableException &e)
         {
             streamlog_out(WARNING) << m_inputTrackerHitRelNames[icol]
                                    << " collection not available" << std::endl;
-            continue;
+            inputHitRels[icol] = nullptr; // optional collection
         }
 
         // reco hit output collections
@@ -193,121 +297,149 @@ void FilterTimeHits::processEvent(LCEvent *evt)
         LCFlagImpl lcFlag(inputHitColls[icol]->getFlag());
         outputTrackerHitColls[icol]->setFlag(lcFlag.getFlag());
 
+        // reco hit contituents output collections, if needed
+        if (hasHitConstituents) {
+            if ((m_outputTrackerHitsConstituentsCollNames[icol] != "") &&
+                (inputHitConstituentsColls[icol] != nullptr))
+            {
+                std::string encoderString = inputHitConstituentsColls[icol]->getParameters().getStringVal("CellIDEncoding");
+                outputTrackerHitConstituentsColls[icol] = new LCCollectionVec(inputHitConstituentsColls[icol]->getTypeName());
+                outputTrackerHitConstituentsColls[icol]->parameters().setValue("CellIDEncoding", encoderString);
+                LCFlagImpl lcFlag(inputHitConstituentsColls[icol]->getFlag());
+                outputTrackerHitConstituentsColls[icol]->setFlag(lcFlag.getFlag());
+            } else {
+                outputTrackerHitConstituentsColls[icol] = nullptr;    
+            }
+        } else {
+            outputTrackerHitConstituentsColls[icol] = nullptr;
+        }
+
         // sim hit output collections
-        outputTrackerSimHitColls[icol] = new LCCollectionVec(inputSimHitColls[icol]->getTypeName());
-        outputTrackerSimHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString);
-        LCFlagImpl lcFlag_sim(inputSimHitColls[icol]->getFlag());
-        outputTrackerSimHitColls[icol]->setFlag(lcFlag_sim.getFlag());
+        if (inputSimHitColls[icol] != nullptr)
+        {
+            outputTrackerSimHitColls[icol] = new LCCollectionVec(inputSimHitColls[icol]->getTypeName());
+            outputTrackerSimHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString);
+            LCFlagImpl lcFlag_sim(inputSimHitColls[icol]->getFlag());
+            outputTrackerSimHitColls[icol]->setFlag(lcFlag_sim.getFlag());
+        }
 
         // reco-sim relation output collections
-        outputTrackerHitRels[icol] = new LCCollectionVec(inputHitRels[icol]->getTypeName());
-        LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag());
-        outputTrackerHitRels[icol]->setFlag(lcFlag_rel.getFlag());
+        if (inputHitRels[icol] != nullptr)
+        {
+            outputTrackerHitRels[icol] = new LCCollectionVec(inputHitRels[icol]->getTypeName());
+            LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag());
+            outputTrackerHitRels[icol]->setFlag(lcFlag_rel.getFlag());
+        }
     }
 
     // --- Loop over the tracker hits and select hits inside the chosen time window:
-    std::vector<std::set<int>> hits_to_save(nTrackerHitCol);
-
+    // IMPORTANT: Cannot assume SimTrkHit and Relation collections contain corresponding elements in the same order. Some sim hits might not have a corrsponding reco hit.
     for (unsigned int icol = 0; icol < inputHitColls.size(); ++icol)
     {
+        // Keep track of (old, new) object pointers for clusters to save in the output collections
+        std::map<TrackerHitPlane *, TrackerHitPlane *> reco_clusters_to_save;
 
         LCCollection *hit_col = inputHitColls[icol];
         if (!hit_col)
+        {
+            streamlog_out(WARNING) << "Cannot retrieve collection: " << m_inputTrackerHitsCollNames[icol] << std::endl;
             continue;
+        }
 
         for (int ihit = 0; ihit < hit_col->getNumberOfElements(); ++ihit)
         {
 
-            if (TrackerHitPlane *hit = dynamic_cast<TrackerHitPlane *>(hit_col->getElementAt(ihit))){
-                // Skipping the hit if its time is outside the acceptance time window
-                double hitT = hit->getTime();
-
-                dd4hep::rec::Vector3D pos = hit->getPosition();
-                double hitR = pos.r();
-
-                // Correcting for the propagation time
-                double dt = hitR / (TMath::C() * m_beta / 1e6);
-                hitT -= dt;
-                streamlog_out(DEBUG3) << "corrected hit at R: " << hitR << " mm by propagation time: " << dt << " ns to T: " << hitT << " ns" << std::endl;
-
-                //Apply time window selection
-                if (hitT < m_time_min || hitT > m_time_max)
-                {
-                    streamlog_out(DEBUG4) << "hit at T: " << hitT << " ns is rejected by timing cuts" << std::endl;
-                    continue;
-                }
-
-                hits_to_save[icol].insert(ihit);
-
-                if (m_fillHistos)
-                    m_corrected_time->Fill(hitT);
+            TrackerHitPlane *hit = dynamic_cast<TrackerHitPlane *>(hit_col->getElementAt(ihit));
+            if (!hit)
+            {
+                streamlog_out(WARNING) << "Cannot retrieve/cast(TrackerHitPlane*) cluster from collection: " << m_inputTrackerHitsCollNames[icol] << std::endl;
+                continue;
             }
+            // Skipping the hit if its time is outside the acceptance time window
+            double hitT = hit->getTime();
 
-        } // ihit loop
+            dd4hep::rec::Vector3D pos = hit->getPosition();
+            double hitR = pos.r();
 
-    } // icol loop
+            // Correcting for the propagation time
+            double dt = hitR / (TMath::C() * m_beta / 1e6);
+            hitT -= dt;
+            streamlog_out(DEBUG3) << "corrected hit at R: " << hitR << " mm by propagation time: " << dt << " ns to T: " << hitT << " ns" << std::endl;
 
-    // --- Add the filtered hits to the output collections:
+            if (m_fillHistos)
+                m_corrected_time_before->Fill(hitT);
 
-    for (unsigned int icol = 0; icol < inputHitColls.size(); ++icol)
-    {
-
-        for (auto &ihit : hits_to_save[icol])
-        {
+            // Apply time window selection
+            if (hitT < m_time_min || hitT > m_time_max)
+            {
+                streamlog_out(DEBUG4) << "hit at T: " << hitT << " ns is rejected by timing cuts" << std::endl;
+                continue;
+            }
 
-            TrackerHitPlane *hit = dynamic_cast<TrackerHitPlane *>(inputHitColls[icol]->getElementAt(ihit));
-            TrackerHitPlaneImpl *hit_new = new TrackerHitPlaneImpl();
-
-            hit_new->setCellID0(hit->getCellID0());
-            hit_new->setCellID1(hit->getCellID1());
-            hit_new->setType(hit->getType());
-            hit_new->setPosition(hit->getPosition());
-            hit_new->setU(hit->getU());
-            hit_new->setV(hit->getV());
-            hit_new->setdU(hit->getdU());
-            hit_new->setdV(hit->getdV());
-            hit_new->setEDep(hit->getEDep());
-            hit_new->setEDepError(hit->getEDepError());
-            hit_new->setTime(hit->getTime());
-            hit_new->setQuality(hit->getQuality());
+            if (m_fillHistos)
+                m_corrected_time_after->Fill(hitT);
 
+            // Save a copy of the cluster into the output collection
+            TrackerHitPlane *hit_new = copyTrackerHitPlane(hit);
             outputTrackerHitColls[icol]->addElement(hit_new);
+            if (outputTrackerHitConstituentsColls[icol])
+            {
+                // save hits contituents as well
+                const EVENT::LCObjectVec &rawHits = hit_new->getRawHits();
+                for (unsigned int ihit = 0; ihit < rawHits.size(); ++ihit)
+                    outputTrackerHitConstituentsColls[icol]->addElement(rawHits[ihit]);
+            }
 
-            LCRelation *rel = dynamic_cast<LCRelation *>(inputHitRels[icol]->getElementAt(ihit));
-
-            SimTrackerHit *simhit = dynamic_cast<SimTrackerHit *>(rel->getTo());
-            SimTrackerHitImpl *simhit_new = new SimTrackerHitImpl();
-
-            simhit_new->setCellID0(simhit->getCellID0());
-            simhit_new->setCellID1(simhit->getCellID1());
-            simhit_new->setPosition(simhit->getPosition());
-            simhit_new->setEDep(simhit->getEDep());
-            simhit_new->setTime(simhit->getTime());
-            simhit_new->setMCParticle(simhit->getMCParticle());
-            simhit_new->setMomentum(simhit->getMomentum());
-            simhit_new->setPathLength(simhit->getPathLength());
-            simhit_new->setQuality(simhit->getQuality());
-            simhit_new->setOverlay(simhit->isOverlay());
-            simhit_new->setProducedBySecondary(simhit->isProducedBySecondary());
-
-            outputTrackerSimHitColls[icol]->addElement(simhit_new);
-
-            LCRelationImpl *rel_new = new LCRelationImpl();
-
-            rel_new->setFrom(hit_new);
-            rel_new->setTo(simhit_new);
-            rel_new->setWeight(rel->getWeight());
-
-            outputTrackerHitRels[icol]->addElement(rel_new);
+            // Keep track of clusters saved
+            reco_clusters_to_save[hit] = hit_new;
 
         } // ihit loop
 
-        streamlog_out(MESSAGE) << " " << hits_to_save[icol].size() << " hits added to the collections: "
-                               << m_outputTrackerHitsCollNames[icol] << ", "
-                               << m_outputTrackerSimHitsCollNames[icol] << ", "
+        // Now loop over the LCRelation collection, and save new SimHits as well as their reco-sim relation
+        LCCollection *hitrel_col = inputHitRels[icol];
+        if (hitrel_col)
+        {
+            for (unsigned int irel = 0; irel < hitrel_col->getNumberOfElements(); ++irel)
+            {
+                LCRelation *rel = dynamic_cast<LCRelation *>(hitrel_col->getElementAt(irel));
+                // now check if the pointer to the reco object (From) is in the list of cluster to save
+                TrackerHitPlane *recohit = dynamic_cast<TrackerHitPlane *>(rel->getFrom());
+                auto reco_hit_itr = reco_clusters_to_save.find(recohit);
+                if (reco_hit_itr != reco_clusters_to_save.end())
+                {
+                    // Reconstructed Hit (already saved in output collection)
+                    TrackerHit *hit_new = reco_hit_itr->second;
+
+                    // Simulated Hit
+                    SimTrackerHit *simhit_new = nullptr;
+                    if (outputTrackerSimHitColls[icol] != nullptr)
+                    {
+                        SimTrackerHit *simhit = dynamic_cast<SimTrackerHit *>(rel->getTo());
+                        simhit_new = copySimTrackerHit(simhit);
+                        outputTrackerSimHitColls[icol]->addElement(simhit_new);
+                    }
+
+                    // LCRelation
+                    LCRelationImpl *rel_new = new LCRelationImpl();
+                    rel_new->setFrom(hit_new);
+                    rel_new->setTo(simhit_new);
+                    rel_new->setWeight(rel->getWeight());
+                    outputTrackerHitRels[icol]->addElement(rel_new);
+                }
+            } // irel loop
+        } // hitrel available
+
+        streamlog_out(MESSAGE) << " " << reco_clusters_to_save.size() << " hits added to the collections: "
+                               << m_outputTrackerHitsCollNames[icol] << ", ";
+        if (m_inputTrackerHitsConstituentsCollNames.size() > 0)
+            if (m_outputTrackerHitsConstituentsCollNames[icol] != "")
+                streamlog_out(MESSAGE) << m_outputTrackerHitsConstituentsCollNames[icol] << ", ";
+        streamlog_out(MESSAGE) << m_outputTrackerSimHitsCollNames[icol] << ", "
                                << m_outputTrackerHitRelNames[icol] << std::endl;
 
         evt->addCollection(outputTrackerHitColls[icol], m_outputTrackerHitsCollNames[icol]);
+        if (outputTrackerHitConstituentsColls[icol])
+            evt->addCollection(outputTrackerHitConstituentsColls[icol], m_outputTrackerHitsConstituentsCollNames[icol]);
         evt->addCollection(outputTrackerSimHitColls[icol], m_outputTrackerSimHitsCollNames[icol]);
         evt->addCollection(outputTrackerHitRels[icol], m_outputTrackerHitRelNames[icol]);
 
@@ -316,8 +448,13 @@ void FilterTimeHits::processEvent(LCEvent *evt)
                               << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type "
                               << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event \n"
                               << " output collection " << m_outputTrackerHitRelNames[icol] << " of type "
-                              << outputTrackerHitRels[icol]->getTypeName() << " added to the event  "
-                              << std::endl;
+                              << outputTrackerHitRels[icol]->getTypeName() << " added to the event  ";
+        if (outputTrackerHitConstituentsColls[icol])
+        {
+            streamlog_out(DEBUG5) << " output collection " << m_outputTrackerHitsConstituentsCollNames[icol] << " of type "
+                                  << outputTrackerHitConstituentsColls[icol]->getTypeName() << " added to the event \n";
+        }
+        streamlog_out(DEBUG5) << std::endl;
 
     } // icol loop