From 3c803ab9038c93d46d3b696da308b1e535029980 Mon Sep 17 00:00:00 2001
From: tmadlener <thomas.madlener@desy.de>
Date: Fri, 16 Feb 2024 15:07:33 +0100
Subject: [PATCH] Make sure to keep necessary type information in output
 relations

---
 .../Digitisers/include/DDTPCDigiProcessor.h   |  9 +++--
 source/Digitisers/src/DDTPCDigiProcessor.cc   | 39 +++++++------------
 2 files changed, 20 insertions(+), 28 deletions(-)

diff --git a/source/Digitisers/include/DDTPCDigiProcessor.h b/source/Digitisers/include/DDTPCDigiProcessor.h
index dfee072..88974a8 100644
--- a/source/Digitisers/include/DDTPCDigiProcessor.h
+++ b/source/Digitisers/include/DDTPCDigiProcessor.h
@@ -58,6 +58,10 @@ class Voxel_tpc;
 
 class TPCModularEndplate ;
 
+namespace UTIL {
+  class LCRelationNavigator;
+}
+
 /** ====== DDTPCDigiProcessor ====== <br>
  *
  * This Processor depends on Circle.h from MarlinUtil
@@ -138,8 +142,8 @@ class DDTPCDigiProcessor : public marlin::Processor {
    */
   virtual void end() ;
   
-  void writeVoxelToHit( Voxel_tpc* aVoxel) ;  
-  void writeMergedVoxelsToHit( std::vector <Voxel_tpc*>* hitList ) ;  
+  void writeVoxelToHit( Voxel_tpc* aVoxel, UTIL::LCRelationNavigator& hitSimHitNav) ;
+  void writeMergedVoxelsToHit( std::vector <Voxel_tpc*>* hitList, UTIL::LCRelationNavigator& hitSimHitNav ) ;
   void plotHelixHitResidual(MCParticle *mcp, CLHEP::Hep3Vector *thisPointRPhi);
   double getPadPhi( CLHEP::Hep3Vector* thisPointRPhi, CLHEP::Hep3Vector* firstPointRPhi, CLHEP::Hep3Vector* middlePointRPhi, CLHEP::Hep3Vector* lastPointRPhi);
   double getPadTheta( CLHEP::Hep3Vector* firstPointRPhi, CLHEP::Hep3Vector* middlePointRPhi, CLHEP::Hep3Vector* lastPointRPhi );
@@ -204,7 +208,6 @@ class DDTPCDigiProcessor : public marlin::Processor {
   int lenpos{};
 
   LCCollectionVec* _trkhitVec{};
-  LCCollectionVec* _relCol{};  
   CellIDEncoder<TrackerHitImpl>* _cellid_encoder {};
 
   int  _NSimTPCHits{};
diff --git a/source/Digitisers/src/DDTPCDigiProcessor.cc b/source/Digitisers/src/DDTPCDigiProcessor.cc
index db91de6..e9157d9 100644
--- a/source/Digitisers/src/DDTPCDigiProcessor.cc
+++ b/source/Digitisers/src/DDTPCDigiProcessor.cc
@@ -32,6 +32,7 @@
 //
 #include "UTIL/LCTrackerConf.h"
 #include <UTIL/ILDConf.h>
+#include <UTIL/LCRelationNavigator.h>
 
 // --- DD4hep ---
 #include "DD4hep/Detector.h"
@@ -500,12 +501,8 @@ void DDTPCDigiProcessor::processEvent( LCEvent * evt )
   
   // created the collection which will be written out 
   _trkhitVec = new LCCollectionVec( LCIO::TRACKERHIT )  ;
-  _relCol = new LCCollectionVec(LCIO::LCRELATION);
-
-  // to store the weights
-  LCFlagImpl lcFlag(0) ;
-  lcFlag.setBit( LCIO::LCREL_WEIGHTED ) ;
-  _relCol->setFlag( lcFlag.getFlag()  ) ;
+  // relations from created trackerhits to the SimTrackerHits that caused them
+  auto hitSimHitNav = UTIL::LCRelationNavigator(LCIO::TRACKERHIT, LCIO::SIMTRACKERHIT);
 
   _cellid_encoder =  new CellIDEncoder<TrackerHitImpl>( lcio::LCTrackerCellID::encoding_string() , _trkhitVec ) ;
   
@@ -1017,7 +1014,7 @@ void DDTPCDigiProcessor::processEvent( LCEvent * evt )
       }
       
       if(seed_hit->getNumberOfAdjacent()==0){ // no adjacent hits so smear and write to hit collection
-        writeVoxelToHit(seed_hit);        
+        writeVoxelToHit(seed_hit, hitSimHitNav);
       }
       
       else if(seed_hit->getNumberOfAdjacent() < (_maxMerge)){ // potential 3-hit cluster, can use simple average merge. 
@@ -1028,7 +1025,7 @@ void DDTPCDigiProcessor::processEvent( LCEvent * evt )
         
         if( clusterSize <= _maxMerge ){ // merge cluster
           seed_hit->setIsMerged();
-          writeMergedVoxelsToHit(hitsToMerge);  
+          writeMergedVoxelsToHit(hitsToMerge, hitSimHitNav);
         }
         delete hitsToMerge;
       } 
@@ -1074,7 +1071,8 @@ void DDTPCDigiProcessor::processEvent( LCEvent * evt )
   
   // add the collection to the event
   evt->addCollection( _trkhitVec , _TPCTrackerHitsCol ) ;
-  evt->addCollection( _relCol , _outRelColName ) ;
+  auto relCol = hitSimHitNav.createLCCollection();
+  evt->addCollection( relCol , _outRelColName ) ;
 
   // delete voxels
   for (unsigned int i = 0; i<_tpcRowHits.size(); ++i){
@@ -1149,7 +1147,7 @@ void DDTPCDigiProcessor::end()
   //  
 }
 
-void DDTPCDigiProcessor::writeVoxelToHit( Voxel_tpc* aVoxel){
+void DDTPCDigiProcessor::writeVoxelToHit( Voxel_tpc* aVoxel, UTIL::LCRelationNavigator& hitSimHitNav){
   
   
   Voxel_tpc* seed_hit  = aVoxel;
@@ -1261,13 +1259,8 @@ void DDTPCDigiProcessor::writeVoxelToHit( Voxel_tpc* aVoxel){
       trkHit->rawHits().push_back( _tpcHitMap[seed_hit] );
     }                        
 
-    LCRelationImpl* rel = new LCRelationImpl;
-    
-    rel->setFrom (trkHit);
-    rel->setTo (_tpcHitMap[seed_hit]);
-    rel->setWeight( 1.0 );
-    _relCol->addElement(rel);
-    
+    hitSimHitNav.addRelation(trkHit, _tpcHitMap[seed_hit], 1.0);
+
     _trkhitVec->addElement( trkHit ); 
     _NRecTPCHits++;
   }
@@ -1298,7 +1291,7 @@ void DDTPCDigiProcessor::writeVoxelToHit( Voxel_tpc* aVoxel){
 #endif
 }
 
-void DDTPCDigiProcessor::writeMergedVoxelsToHit( vector <Voxel_tpc*>* hitsToMerge){
+void DDTPCDigiProcessor::writeMergedVoxelsToHit( vector <Voxel_tpc*>* hitsToMerge, UTIL::LCRelationNavigator& hitSimHitNav){
   
   
   TrackerHitImpl* trkHit = new TrackerHitImpl ;
@@ -1324,13 +1317,9 @@ void DDTPCDigiProcessor::writeMergedVoxelsToHit( vector <Voxel_tpc*>* hitsToMerg
       trkHit->rawHits().push_back( _tpcHitMap[hitsToMerge->at(ihitCluster)] );
     }                        
 
-    LCRelationImpl* rel = new LCRelationImpl;
-    
-    rel->setFrom (trkHit);
-    rel->setTo (_tpcHitMap[ hitsToMerge->at(ihitCluster) ]);
-    rel->setWeight( float(1.0/number_of_hits_to_merge) );
-    _relCol->addElement(rel);
-    
+    hitSimHitNav.addRelation(trkHit,
+                             _tpcHitMap[hitsToMerge->at(ihitCluster)],
+                             float(1.0/number_of_hits_to_merge));
   }
   
   double avgZ = sumZ/(hitsToMerge->size());