From 6c695987ae0dae18fd52239020c419f025298b4b Mon Sep 17 00:00:00 2001
From: Chengdong Fu <fucd@ihep.ac.cn>
Date: Fri, 23 Oct 2020 18:38:05 +0800
Subject: [PATCH] update and fix runtime error on CEPCSW sim

---
 .../Clupatra => Tracking}/TrackingHelper.h    |   6 +-
 .../Tracking/src/Clupatra/ClupatraAlg.cpp     | 190 ++++++++++--------
 .../Tracking/src/Clupatra/ClupatraAlg.h       |  11 +-
 .../Tracking/src/Clupatra/clupatra_new.cpp    |  43 ++--
 .../Tracking/src/Clupatra/clupatra_new.h      |  28 +--
 5 files changed, 146 insertions(+), 132 deletions(-)
 rename Reconstruction/Tracking/{src/Clupatra => Tracking}/TrackingHelper.h (92%)

diff --git a/Reconstruction/Tracking/src/Clupatra/TrackingHelper.h b/Reconstruction/Tracking/Tracking/TrackingHelper.h
similarity index 92%
rename from Reconstruction/Tracking/src/Clupatra/TrackingHelper.h
rename to Reconstruction/Tracking/Tracking/TrackingHelper.h
index 2a16aa82..d3bcd2ef 100644
--- a/Reconstruction/Tracking/src/Clupatra/TrackingHelper.h
+++ b/Reconstruction/Tracking/Tracking/TrackingHelper.h
@@ -7,6 +7,7 @@
 #include "UTIL/CellIDDecoder.h"
 #include "UTIL/ILDConf.h"
 #include "lcio.h"
+#include <array>
 
 inline bool hasTrackStateAt(edm4hep::ConstTrack track, int location) {
     for (auto it = track.trackStates_begin(); it != track.trackStates_end(); it++) {
@@ -20,12 +21,15 @@ inline bool hasTrackStateAt(edm4hep::ConstTrack track, int location) {
 inline edm4hep::TrackState getTrackStateAt(edm4hep::ConstTrack track, int location) {
     for (auto it = track.trackStates_begin(); it != track.trackStates_end(); it++) {
         if (it->location == location) {
-            return *it;
+	  return *it;
         }
     }
     return edm4hep::TrackState();
 }
 
+inline std::array<float,15> getCovMatrix(const edm4hep::ConstTrack &track) {
+    return track.getTrackStates(0).covMatrix;
+}
 inline float getTanLambda(const edm4hep::ConstTrack &track) {
     return track.getTrackStates(0).tanLambda;
 }
diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
index c048c098..af73dccf 100644
--- a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
+++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
@@ -15,13 +15,13 @@
 //---- MarlinUtil
 
 //---- LCIO ---
-#include "IMPL/LCCollectionVec.h"
-#include "EVENT/SimTrackerHit.h"
-#include "UTIL/Operators.h"
-#include "UTIL/LCTOOLS.h"
-#include "UTIL/CellIDDecoder.h"
+//#include "IMPL/LCCollectionVec.h"
+//#include "EVENT/SimTrackerHit.h"
+//#include "UTIL/Operators.h"
+//#include "UTIL/LCTOOLS.h"
+//#include "UTIL/CellIDDecoder.h"
 #include "UTIL/ILDConf.h"
-#include "UTIL/LCIterator.h"
+//#include "UTIL/LCIterator.h"
 
 //-------gsl -----
 #include "gsl/gsl_randist.h"
@@ -198,29 +198,32 @@ void ClupatraAlg::printParameters() {
 StatusCode ClupatraAlg::initialize() {
 
 	// Usually a good idea to
-	printParameters() ;
+        // don't need, since Gaudi Algorithm will print all Property  
+	//printParameters() ;
 
 	auto _trackSystemSvc = service<ITrackSystemSvc>("TrackSystemSvc");
 	if ( !_trackSystemSvc ) {
 		error() << "Fail to find TrackSystemSvc ..." << endmsg;
 	}
-	// _trksystem =  MarlinTrk::Factory::createMarlinTrkSystem( "KalTest" , marlin::Global::GEAR , "" ) ;
-
-	_trksystem = _trackSystemSvc->getTrackSystem();
-	// FIXME Mingrui Important fix it later
-	// _trksystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::useQMS,        true ) ;
-	_trksystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::useQMS,        false ) ;
-	_trksystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::usedEdx,       true) ;
-	_trksystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::useSmoothing,  false) ;
+
+	_trksystem = _trackSystemSvc->getTrackSystem(this);
+	if(!_trksystem){
+	  error() << "Cannot initialize MarlinTrkSystem of Type: KalTest" <<endmsg;
+	  return StatusCode::FAILURE;
+	}
+
+	_trksystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::useQMS,        _MSOn ) ;
+	_trksystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::usedEdx,       _ElossOn) ;
+	_trksystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::useSmoothing,  _SmoothOn) ;
 	_trksystem->init() ;
 
 	_nRun = 0 ;
 	_nEvt = 0 ;
-
-        tree = new TTree("Tuple", "Particle Tree");
-        tree->Branch("omega", &omega, "omega/D");
-        tree->Branch("totalCandidates", &totalCandidates, "totalCandidates/I");
-        tree->Branch("eventNumber", &_nEvt, "eventNumber/I");
+	// FIXME: fucd
+        //tree = new TTree("Tuple", "Particle Tree");
+        //tree->Branch("omega", &omega, "omega/D");
+        //tree->Branch("totalCandidates", &totalCandidates, "totalCandidates/I");
+        //tree->Branch("eventNumber", &_nEvt, "eventNumber/I");
 
 	return GaudiAlgorithm::initialize();
 
@@ -299,7 +302,7 @@ StatusCode ClupatraAlg::execute() {
 	//   create clupa and clustering hits for every edm4hep hit
 	//===============================================================================================
 	if (col == nullptr) {
-                return StatusCode::SUCCESS;
+	  return StatusCode::SUCCESS;
         }
 
 
@@ -308,8 +311,6 @@ StatusCode ClupatraAlg::execute() {
 	clupaHits.resize( nHit ) ;       // creates clupa hits (w/ default c'tor)
 	nncluHits.reserve( nHit ) ;
 
-
-	// FIXME: Mingrui fix the output message
 	debug()  << "  create clupatra TPC hits, n = " << nHit << endmsg ;
 
 	// Mingrui!: Copy the items in col to col_copy to avoid losing the pointer.
@@ -317,7 +318,7 @@ StatusCode ClupatraAlg::execute() {
 
 		//------
 		ConstTrackerHit th(col->at(i));
-		// debug() << i << " " << th->getCellID0() << " " << th->getCellID1() << endmsg;
+		//debug() << i << " " << th->getCellID() << endmsg;
 		if ( fabs(th.getPosition()[2]) > driftLength ) continue;
 
 		ClupaHit* ch  = & clupaHits[i] ;
@@ -339,8 +340,7 @@ StatusCode ClupatraAlg::execute() {
 		// May cause a problem
 		ch->layer = getLayer( th );
 
-		// FIXME: Mingrui Debug
-		// streamlog_out( DEBUG ) << "  ch->layer = idDec( th )[ ILDCellID0::layer ] = " <<  ch->layer << " - CellID0 " << th->getCellID0() << std::endl ;
+		//debug() << "ch->layer = idDec( th )[ ILDCellID0::layer ] = " <<  ch->layer << " - CellID0 " << th.getCellID() << endmsg;
 
 		ch->zIndex = zIndex( &th ) ;
 
@@ -358,7 +358,7 @@ StatusCode ClupatraAlg::execute() {
 	HitListVector hitsInLayer( maxTPCLayers ) ;
 	addToHitListVector(  nncluHits.begin(), nncluHits.end() , hitsInLayer  ) ;
 
-	// debug() << "  added  " <<  nncluHits.size()  << "  tp hitsInLayer - > size " <<  hitsInLayer.size() << std::endl ;
+	debug() << "  added  " <<  nncluHits.size()  << "  tp hitsInLayer - > size " <<  hitsInLayer.size() << endmsg;
 
 	//---------------------------------------------------------------------------------------------------------
 
@@ -403,9 +403,9 @@ StatusCode ClupatraAlg::execute() {
 
 	TrackCollection* outCol =  _ClupatraTrackCollectionHandle.createAndPut();
 	TrackCollection* tsCol  =  _ClupatraTrackSegmentCollectionHandle.createAndPut();
-        std::vector<ClupaPlcioConstTrack*> outCol_tmp;
-        std::vector<ClupaPlcioConstTrack*> tsCol_tmp;
-
+        std::vector<ClupaPlcioTrack*> outCol_tmp;
+        std::vector<ClupaPlcioTrack*> tsCol_tmp;
+	
 	//---------------------------------------------------------------------------------------------------------
 
 	timer.time(t_init ) ;
@@ -429,9 +429,9 @@ StatusCode ClupatraAlg::execute() {
 	IMarlinTrkFitter fitter( _trksystem ) ;
 
 
-	debug() << "===============================================================================================\n"
-		<< "   first step of Clupatra algorithm: find seeds with NN-clustering  in " <<  _nLoop << " loops - max dist = " << _distCut <<" \n"
-		<< "===============================================================================================\n"  << endmsg;
+	debug() << "===============================================================================================" << endmsg;
+	debug() << "   first step of Clupatra algorithm: find seeds with NN-clustering  in " <<  _nLoop << " loops - max dist = " << _distCut << endmsg;
+	debug() << "==============================================================================================="  << endmsg;
 
 	// ---- introduce a loop over increasing distance cuts for finding the tracks seeds
 	//      -> should fix (some of) the problems seen @ 3 TeV with extremely boosted jets
@@ -819,8 +819,8 @@ StatusCode ClupatraAlg::execute() {
 
 		MarlinTrk::IMarlinTrack* trk = fit( *icv ) ;
 		trk->smooth() ;
-		edm4hep::ConstTrack edm4hepTrk = converter( *icv ) ;
-		tsCol_tmp.push_back( new ClupaPlcioConstTrack(edm4hepTrk) ) ;
+		edm4hep::Track edm4hepTrk = converter( *icv ) ;
+		tsCol_tmp.push_back( new ClupaPlcioTrack(edm4hepTrk) ) ;
 		MarTrk_of_edm4hepTrack(edm4hepTrk) = 0 ;
 		delete trk ;
 	}
@@ -836,14 +836,13 @@ StatusCode ClupatraAlg::execute() {
 	//  compute some track parameters for possible merging
 	//===============================================================================================
 
-	typedef nnclu::NNClusterer<ClupaPlcioConstTrack> TrackClusterer ;
+	typedef nnclu::NNClusterer<ClupaPlcioTrack> TrackClusterer ;
 	TrackClusterer nntrkclu ;
-	MakePLCIOElement<ClupaPlcioConstTrack> trkMakeElement ;
+	MakePLCIOElement<ClupaPlcioTrack> trkMakeElement ;
 
 	for( int i=0,N=tsCol_tmp.size() ;  i<N ; ++i ) {
-
-		ConstTrack track = tsCol_tmp.at(i)->edm4hepTrack;
-		computeTrackInfo(track) ;
+	  edm4hep::ConstTrack track = tsCol_tmp.at(i)->edm4hepTrack;
+	  computeTrackInfo(track) ;
 	}
 
 	//===============================================================================================
@@ -867,7 +866,7 @@ StatusCode ClupatraAlg::execute() {
 
 			for( int i=0,N=tsCol_tmp.size() ;  i<N ; ++i ){
 
-				edm4hep::ConstTrack trk = tsCol_tmp.at(i)->edm4hepTrack;
+			        edm4hep::ConstTrack trk = tsCol_tmp.at(i)->edm4hepTrack;
 
 				const TrackInfoStruct* ti = TrackInfo_of_edm4hepTrack(trk);
 
@@ -916,7 +915,7 @@ StatusCode ClupatraAlg::execute() {
 
 					//streamlog_out( DEBUG3 ) << lcshort(  (*itC)->first ) << std::endl ;
 
-					edm4hep::ConstTrack trk = (*itC)->first->edm4hepTrack ;
+					edm4hep::ConstTrack trk = (*itC)->first->edm4hepTrack;
 
 					mergedTrk.push_back( trk ) ;
 
@@ -953,8 +952,8 @@ StatusCode ClupatraAlg::execute() {
 
 				MarlinTrk::IMarlinTrack* mTrk = fit( &hits ) ;
 				mTrk->smooth() ;
-				edm4hep::ConstTrack track = converter( &hits ) ;
-				tsCol_tmp.push_back( new ClupaPlcioConstTrack(track) ) ;
+				edm4hep::Track track = converter( &hits ) ;
+				tsCol_tmp.push_back( new ClupaPlcioTrack(track) ) ;
 				MarTrk_of_edm4hepTrack(track) = 0 ;
 				delete mTrk ;
 				computeTrackInfo( track ) ;
@@ -982,10 +981,10 @@ StatusCode ClupatraAlg::execute() {
 		TrackClusterer::cluster_vector curSegCluVec ;
 		curSegCluVec.setOwner() ;
 
-
+		//for(std::vector<ClupaPlcioTrack*>::iterator it=tsCol_tmp.begin();it!=tsCol_tmp.end();it++){
 		for( int i=tsCol_tmp.size()-1 ;  i>=0 ; --i ){
 
-			ConstTrack trk = tsCol_tmp.at(i)->edm4hepTrack;
+		        edm4hep::Track trk = tsCol_tmp.at(i)->edm4hepTrack;
 
 			std::bitset<32> type = trk.getType() ;
 
@@ -1006,18 +1005,21 @@ StatusCode ClupatraAlg::execute() {
 			} else {   // ... is not a curler ->  add a copy to the final tracks collection
 
 
-				if( copyTrackSegments) {
-
-					outCol_tmp.push_back( new ClupaPlcioConstTrack(trk)  ) ;
+			        if( copyTrackSegments) {
+				  
+				        outCol_tmp.push_back( new ClupaPlcioTrack(trk)  ) ;
 
 				}else{
 
-					outCol_tmp.push_back( new ClupaPlcioConstTrack(trk) ) ;
-
+				        outCol_tmp.push_back( new ClupaPlcioTrack(trk) ) ;
+					debug() << "now outCol_tmp size = " << outCol_tmp.size() << endmsg;
 					// FIXME !!! Possibly different with the LCIO version
 					// tsCol-> Do not remove it at the Beginning
 					// Very very important
 					// tsCol->removeElementAt( i ) ;
+					// fucd
+					//delete tsCol_tmp.at(i);
+					tsCol_tmp.erase(tsCol_tmp.begin()+i);
 				}
 
 				// FIXME Mingrui remove
@@ -1032,12 +1034,11 @@ StatusCode ClupatraAlg::execute() {
 
 		nntrkclu.cluster( curSegVec.begin() , curSegVec.end() , std::back_inserter( curSegCluVec ), trkMerge , 2  ) ;
 
-		// FIXME: Mingrui
 		debug() << " ===== merged tracks - # cluster: " << curSegCluVec.size()
 			<< " from " << tsCol_tmp.size() << " track segments "    << "  ============================== " << endmsg ;
 
-for (auto t : tsCol_tmp)
-// std::cout << t->edm4hepTrack << std::endl;
+		//for (auto t : tsCol_tmp)
+		//debug() << t->edm4hepTrack << endmsg;
 
 // debug() << "Finish this line" << endmsg;
 
@@ -1048,14 +1049,13 @@ for (auto t : tsCol_tmp)
 
 			TrackClusterer::cluster_type*  curSegClu = *it ;
 
-			std::list<ConstTrack> mergedTrk ;
+			std::list<edm4hep::Track> mergedTrk ;
 
 			for( TrackClusterer::cluster_type::iterator itC = curSegClu->begin() ; itC != curSegClu->end() ; ++ itC ){
 
-				// FIXME: Mingrui
-				// streamlog_out( DEBUG4 ) << lcshort(  (*itC)->first ) << std::endl ;
-
-				mergedTrk.push_back( (*itC)->first->edm4hepTrack ) ;
+			  //debug() << lcshort(  (*itC)->first ) << endmsg;
+			  //debug() << getOmega((*itC)->first->edm4hepTrack) << endmsg;
+			  mergedTrk.push_back( (*itC)->first->edm4hepTrack ) ;
 			}
 
 
@@ -1074,7 +1074,7 @@ for (auto t : tsCol_tmp)
 
 				// == and copy all the hits
 				unsigned hitCount = 0 ;
-				for( std::list<ConstTrack>::iterator itML = mergedTrk.begin() ; itML != mergedTrk.end() ; ++ itML ){
+				for( std::list<edm4hep::Track>::iterator itML = mergedTrk.begin() ; itML != mergedTrk.end() ; ++ itML ){
 
 					for (auto itHit = (*itML).trackerHits_begin(); itHit != (*itML).trackerHits_end(); itHit++) {
 						trk.addToTrackerHits(*itHit);
@@ -1141,36 +1141,39 @@ for (auto t : tsCol_tmp)
 
 
 
-				// streamlog_out( DEBUG2 ) << "   create new merged track from bestTrack parameters - ptr to MarlinTrk : " << trk->ext<MarTrk>()
-				//  << "   with subdetector hit numbers  : " <<  trk->subdetectorHitNumbers()[0 ] << " , " <<  trk->subdetectorHitNumbers()[1]
-				//  << std::endl ;
+				//debug() << "   create new merged track from bestTrack parameters - ptr to MarlinTrk : " << trk->ext<MarTrk>()
+				//	<< "   with subdetector hit numbers  : " <<  trk->subdetectorHitNumbers()[0 ] << " , " <<  trk->subdetectorHitNumbers()[1]
+				//	<< endmsg;
 
 
-				outCol_tmp.push_back( new ClupaPlcioConstTrack(trk) )  ;
+				outCol_tmp.push_back( new ClupaPlcioTrack(trk) )  ;
 
 			} else { //==========================
 
 				// we move the first segment to the final list and keep pointers to the other segments
 
-				std::list<ConstTrack>::iterator itML = mergedTrk.begin() ;
+			        std::list<edm4hep::Track>::iterator itML = mergedTrk.begin() ;
 
-				ConstTrack trk = (ConstTrack) *itML++ ;
+				edm4hep::Track trk = (edm4hep::Track) *itML++ ;
 
 				for(  ; itML != mergedTrk.end() ; ++itML ){
 
 					// add a pointer to the original track segment
 					// FIXME I need to add sub track
-					// trk->addTrack(**itML) ;
+				  trk.addToTracks(*itML) ;
 				}
 
-				outCol_tmp.push_back( new ClupaPlcioConstTrack(trk) ) ;
-
+				outCol_tmp.push_back( new ClupaPlcioTrack(trk) ) ;
+				debug() << "now outCol_tmp size = " << outCol_tmp.size() << endmsg;
 				//remove from segment collection:
 				for( int i=0,N=tsCol_tmp.size() ; i<N ; ++i) {
-					if( ( tsCol_tmp.at(i)->edm4hepTrack ) == trk ){
-						// FIXME Do not remove
-						// tsCol_tmp.removeElementAt( i ) ;
-						break ;
+				  if( ( tsCol_tmp.at(i)->edm4hepTrack ) == trk ){
+					        // FIXME Do not remove
+						//tsCol_tmp.removeElementAt( i ) ;
+					  //fucd
+					  //delete tsCol_tmp[i];
+					  tsCol_tmp.erase(tsCol_tmp.begin()+i);
+					  break ;
 					}
 				}
 
@@ -1180,7 +1183,7 @@ for (auto t : tsCol_tmp)
 		//---------------------------------------------------------------------------------------------
 		// // add all tracks that have not been merged :
 
-// debug() << "Still works here" << endmsg;
+		debug() << "Still works here" << endmsg;
 		for( TrackClusterer::element_vector::iterator it = curSegVec.begin(); it != curSegVec.end() ;++it){
 
                                 // std::cout << "Can I get the track? 1" << std::endl;
@@ -1188,13 +1191,13 @@ for (auto t : tsCol_tmp)
 
                                 // std::cout << "Can I get the track? 2" << std::endl;
                                 // std::cout << (*it)->first->edm4hepTrack;
-				ConstTrack trk = (*it)->first->edm4hepTrack ;
+			        edm4hep::Track trk = (*it)->first->edm4hepTrack ;
                                 // std::cout << "Can I get the track? 3" << std::endl;
                                 // std::cout << trk << std::endl;
 
 				if( copyTrackSegments) {
 
-					ConstTrack t =  trk;
+				        edm4hep::Track t =  trk;
 
 					TrackInfo_of_edm4hepTrack(t) = 0 ; // set extension to 0 to prevent double free ...
 
@@ -1203,17 +1206,20 @@ for (auto t : tsCol_tmp)
 					// FIXME Mingrui debug
 					// streamlog_out( DEBUG2 ) << "   create new track from existing LCIO track  - ptr to MarlinTrk : " << t->ext<MarTrk>()  << std::endl ;
 
-					outCol_tmp.push_back( new ClupaPlcioConstTrack(t) ) ;
+					outCol_tmp.push_back( new ClupaPlcioTrack(t) ) ;
 
 				} else {
-					outCol_tmp.push_back( new ClupaPlcioConstTrack(trk) ) ;
-
+				        outCol_tmp.push_back( new ClupaPlcioTrack(trk) ) ;
+				        debug() << "now outCol_tmp size = " << outCol_tmp.size() << endmsg;
 					//remove from segment collection:
 					for( int i=0,N=tsCol_tmp.size() ; i<N ; ++i) {
-						if( ( tsCol_tmp.at(i)->edm4hepTrack ) == trk ){
-							// FIXME Do not remove
+					  if( ( tsCol_tmp.at(i)->edm4hepTrack ) == trk ){
+						        // FIXME Do not remove
 							// tsCol->removeElementAt( i ) ;
-							break ;
+						  //fucd
+						  //delete tsCol_tmp[i];
+						  tsCol_tmp.erase(tsCol_tmp.begin()+i);
+						  break ;
 						}
 					}
 				}
@@ -1241,13 +1247,18 @@ for (auto t : tsCol_tmp)
 	timer.time( t_pickup ) ;
 
         totalCandidates = outCol_tmp.size();
+        debug() << "Total " << totalCandidates << " clupatra tracks to save" << endmsg;
         for (auto t : outCol_tmp) {
-            // std::cout << "Final track: omega = " << getOmega(t->edm4hepTrack) << std::endl;
-            // omega = getOmega(t->edm4hepTrack);
-            // tree->Fill();
-            outCol->push_back(t->edm4hepTrack);
+	  debug() << "Final track: omega = " << getOmega(t->edm4hepTrack) << endmsg;
+	  // omega = getOmega(t->edm4hepTrack);
+	  // tree->Fill();
+	  outCol->push_back(t->edm4hepTrack);
         }
-        for (auto t : outCol_tmp) {
+	for (auto t : tsCol_tmp) {
+	  debug() << "Final track segment: omega = " << getOmega(t->edm4hepTrack) << endmsg;
+	  tsCol->push_back(t->edm4hepTrack);
+	}
+	for (auto t : outCol_tmp) {
             delete t;
         }
         for (auto t : tsCol_tmp) {
@@ -1314,7 +1325,8 @@ void ClupatraAlg::computeTrackInfo(  edm4hep::ConstTrack lTrk  ){
 	ti->zMin = zMin ;
 	ti->zMax = zMax ;
 	ti->zAvg = zAvg ;
-
+	
+	//debug() << lTrk.id() << " " << tsF.D0 << " " << tsF.phi << " " << tsF.omega << " " << tsF.Z0 << " " << tsF.tanLambda << endmsg;
 }
 
 
@@ -1361,7 +1373,7 @@ streamlog_out( DEBUG5 ) <<   " ------------------- check()  done " << std::endl
 //====================================================================================================
 
 StatusCode ClupatraAlg::finalize(){
-        tree->SaveAs("Tuple.root");
+  //tree->SaveAs("Tuple.root");
 
 	/*
 	 * FIXME debug
diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
index 80c10b64..7cdbb0f4 100644
--- a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
+++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
@@ -15,7 +15,7 @@
 #include "gear/TPCModule.h"
 #include "RuntimeMap.h"
 #include "clupatra_new.h"
-#include "TrackingHelper.h"
+#include "Tracking/TrackingHelper.h"
 
 #include "TrackSystemSvc/MarlinTrkUtils.h"
 #include "TrackSystemSvc/HelixTrack.h"
@@ -135,7 +135,9 @@ class ClupatraAlg : public GaudiAlgorithm {
   Gaudi::Property<float> _trackEndsOuterCentralDist{this, "TrackEndsOuterCentralDist", 25.}; // "maximum radial distance [mm] from outer field cage of last hit, such that the track is considered to end at the end " 
   Gaudi::Property<float> _trackEndsOuterForwardDist{this, "TrackEndsOuterForwardDist" , 40.}; // "maximum distance in z [mm] from endplate of last hit, such that the track is considered to end at the end " ,
   Gaudi::Property<float> _trackIsCurlerOmega{this, "TrackIsCurlerOmega", 0.001}; // "minimum curvature omega of a track segment for being considered a curler"
-
+  Gaudi::Property<bool> _MSOn{this, "MultipleScatteringOn", false};
+  Gaudi::Property<bool> _ElossOn{this, "EnergyLossOn", true};
+  Gaudi::Property<bool> _SmoothOn{this, "SmoothOn", false};
 
 
   DataHandle<edm4hep::TrackerHitCollection> _TPCHitCollectionHandle{"TPCTrackerHits", Gaudi::DataHandle::Reader, this};
@@ -147,11 +149,6 @@ class ClupatraAlg : public GaudiAlgorithm {
 
   float _bfield ;
 
-  bool _MSOn ;
-  bool _ElossOn ;
-  bool _SmoothOn ;
-
-
   int _nRun ;
   int _nEvt ;
 
diff --git a/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp b/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
index 479b981c..0751a4d0 100644
--- a/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
+++ b/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
@@ -29,18 +29,18 @@ namespace lcio{
 
 extern gear::GearMgr* gearMgr; // = _gear->getGearMgr();
 extern RuntimeMap<clupatra_new::CluTrack*, MarlinTrk::IMarlinTrack*> MarTrkof;
-extern RuntimeMap<edm4hep::ConstTrack, MarlinTrk::IMarlinTrack*> MarTrk_of_edm4hepTrack;
-extern RuntimeMap<edm4hep::ConstTrack, clupatra_new::TrackInfoStruct*> TrackInfo_of_edm4hepTrack;
+extern RuntimeMap<edm4hep::Track, MarlinTrk::IMarlinTrack*> MarTrk_of_edm4hepTrack;
+extern RuntimeMap<edm4hep::Track, clupatra_new::TrackInfoStruct*> TrackInfo_of_edm4hepTrack;
 
 namespace clupatra_new{
 
-	bool TrackZSort::operator()( edm4hep::ConstTrack l, edm4hep::ConstTrack r){
+	bool TrackZSort::operator()( edm4hep::Track l, edm4hep::Track r){
 		return (  std::abs( TrackInfo_of_edm4hepTrack(l)->zAvg )   <   std::abs( TrackInfo_of_edm4hepTrack(r)->zAvg )  ) ;
 	}
 
-	void ComputeTrackerInfo::operator()( edm4hep::ConstTrack o )
+	void ComputeTrackerInfo::operator()( edm4hep::Track o )
 	{
-		edm4hep::ConstTrack lTrk  = o ;
+		edm4hep::Track lTrk  = o ;
 
 		// compute z-extend of this track segment
 		// const edm4hep::TrackerHitVec& hv = lTrk->getTrackerHits() ;
@@ -69,10 +69,10 @@ namespace clupatra_new{
 		TrackInfo_of_edm4hepTrack(lTrk)->zAvg = zAvg ;
 
 	}
-	bool TrackCircleDistance:: operator()( nnclu::Element<ClupaPlcioConstTrack>* h0, nnclu::Element<ClupaPlcioConstTrack>* h1){
+	bool TrackCircleDistance:: operator()( nnclu::Element<ClupaPlcioTrack>* h0, nnclu::Element<ClupaPlcioTrack>* h1){
 
-		edm4hep::ConstTrack trk0 = h0->first->edm4hepTrack ;
-		edm4hep::ConstTrack trk1 = h1->first->edm4hepTrack ;
+		edm4hep::Track trk0 = h0->first->edm4hepTrack ;
+		edm4hep::Track trk1 = h1->first->edm4hepTrack ;
 
 		const TrackInfoStruct* ti0 =  TrackInfo_of_edm4hepTrack(trk0);
 		const TrackInfoStruct* ti1 =  TrackInfo_of_edm4hepTrack(trk1);
@@ -1300,8 +1300,8 @@ start:
 
 	//---------------------------------------------------------------------------------------------------------------------------
 
-	edm4hep::ConstTrack PLCIOTrackConverter::operator() (CluTrack* c) {
-
+	edm4hep::Track PLCIOTrackConverter::operator() (CluTrack* c) {
+	  
 		static lcio::BitField64 encoder( lcio::ILDCellID0::encoder_string ) ;
 
 		edm4hep::Track trk;
@@ -1315,8 +1315,8 @@ start:
 			// reset outliers (not used in fit)  bit
 			//      IMPL::TrackerHitImpl* thi = dynamic_cast<IMPL::TrackerHitImpl*> (  (*hi)->first->edm4hepHit )  ;
 			//      thi->setQualityBit( UTIL::ILDTrkHitQualityBit::USED_IN_FIT , 0 )  ;
-
-			trk.addToTrackerHits(  (*hi)->first->edm4hepHit ) ;
+		  
+		        trk.addToTrackerHits(  (*hi)->first->edm4hepHit ) ;
 			e += (*hi)->first->edm4hepHit.getEDep() ;
 			nHit++ ;
 		}
@@ -1379,8 +1379,6 @@ start:
 				edm4hep::ConstTrackerHit fHit = (hitsInFit.back().first);
 				edm4hep::ConstTrackerHit lHit = (hitsInFit.front().first);
 
-
-
 				//order of hits in fit is reversed wrt time  (we fit inwards)
 
 				// ======= get TrackState at first hit  ========================
@@ -1402,11 +1400,16 @@ start:
 				code = mtrk->getTrackState( lHit, tsLH, chi2, ndf ) ;
 #else     // get the track state at the last hit by propagating from the last(first) constrained fit position (a la MarlinTrkUtils)
 				edm4hep::ConstTrackerHit last_constrained_hit;
-				mtrk->getTrackerHitAtPositiveNDF( last_constrained_hit );
+				code = mtrk->getTrackerHitAtPositiveNDF( last_constrained_hit );
+				mtrk->smooth() ;
+				if( code != MarlinTrk::IMarlinTrack::success ){
+				  std::cout << "last_constrained_hit is unavaibale" << std::endl;
+				}
+				else{
 // std::cout << "the hit" << std::endl;
 // std::cout << lHit.getCellID0() << std::endl;
 // std::cout << last_constrained_hit << std::endl;
-				code = mtrk->smooth() ;
+				//code = mtrk->smooth() ;
 // std::cout << "Smooth success" << std::endl;
 				// gear::Vector3D last_hit_pos( lHit->getPosition()[0], lHit->getPosition()[1], lHit->getPosition()[2] );
 // std::cout << "lHit = " << lHit << std::endl;
@@ -1415,13 +1418,11 @@ start:
 // std::cout << "Position" << lHit.getPosition() << std::endl;
 // std::cout << "----------------------------------------------------------" << std::endl;
 
-
-				edm4hep::Vector3d last_hit_pos( lHit.getPosition() );
+				  edm4hep::Vector3d last_hit_pos( lHit.getPosition() );
 // std::cout << "lHit = " << lHit << std::endl;
-				code = mtrk->propagate( last_hit_pos, last_constrained_hit, tsLH, chi2, ndf);
-
+				  code = mtrk->propagate( last_hit_pos, last_constrained_hit, tsLH, chi2, ndf);
+                                } 
 // std::cout << "Propagate success" << std::endl;
-
 #endif
 
 // std::cout << "We are here" << std::endl; 
diff --git a/Reconstruction/Tracking/src/Clupatra/clupatra_new.h b/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
index 5e5c9529..e8aecc0c 100644
--- a/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
+++ b/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
@@ -2,7 +2,7 @@
 #define clupatra_new_h
 
 #include "RuntimeMap.h"
-#include "TrackingHelper.h"
+#include "Tracking/TrackingHelper.h"
 
 #include <cmath>
 #include <algorithm>
@@ -81,9 +81,9 @@ namespace clupatra_new{
 	typedef std::list<Hit*>        HitList ;
 	typedef std::vector< HitList > HitListVector ;
 
-	struct ClupaPlcioConstTrack {
-		edm4hep::ConstTrack edm4hepTrack ;
-                ClupaPlcioConstTrack(edm4hep::ConstTrack edm4hepTrack) : edm4hepTrack(edm4hepTrack) {}
+	struct ClupaPlcioTrack {
+		edm4hep::Track edm4hepTrack ;
+                ClupaPlcioTrack(edm4hep::Track edm4hepTrack) : edm4hepTrack(edm4hepTrack) {}
 	};
 
 	// typedef GenericHitVec<ClupaHit>      GHitVec ;
@@ -138,7 +138,7 @@ namespace clupatra_new{
 	//------------------------------------------------------------------------------------------
 
 	struct PtSort {  // sort tracks wtr to pt - largest first
-		inline bool operator()( const edm4hep::ConstTrack l, const edm4hep::ConstTrack r) {
+		inline bool operator()( const edm4hep::Track l, const edm4hep::Track r) {
 			return ( std::abs( getOmega(l) ) < std::abs( getOmega(r) )  );  // pt ~ 1./omega
 		}
 	};
@@ -222,7 +222,7 @@ namespace clupatra_new{
 		bool UsePropagate ;
 		PLCIOTrackConverter() : UsePropagate(false ) {}
 
-		edm4hep::ConstTrack operator() (CluTrack* c) ;
+		edm4hep::Track operator() (CluTrack* c) ;
 
 	} ;
 
@@ -354,7 +354,7 @@ namespace clupatra_new{
 	/** Helper class to compute track segment properties.
 	*/
 	struct ComputeTrackerInfo{
-		void operator()( edm4hep::ConstTrack o );
+		void operator()( edm4hep::Track o );
 	};
 
 	//=======================================================================================
@@ -372,10 +372,10 @@ namespace clupatra_new{
 			float _b ;
 
 			/** Merge condition: ... */
-			inline bool operator()( nnclu::Element<ClupaPlcioConstTrack>* h0, nnclu::Element<ClupaPlcioConstTrack>* h1){
+			inline bool operator()( nnclu::Element<ClupaPlcioTrack>* h0, nnclu::Element<ClupaPlcioTrack>* h1){
 
-				edm4hep::ConstTrack trk0 = h0->first->edm4hepTrack ;
-				edm4hep::ConstTrack trk1 = h1->first->edm4hepTrack ;
+				edm4hep::Track trk0 = h0->first->edm4hepTrack ;
+				edm4hep::Track trk1 = h1->first->edm4hepTrack ;
 
 
 				// protect against merging multiple segments (and thus complete tracks)
@@ -423,8 +423,8 @@ namespace clupatra_new{
 
 				// now we take the larger segment and see if we can add the three hits from the other segment...
 
-				edm4hep::ConstTrack trk = ( nhit0 > nhit1 ? trk0 :  trk1 ) ;
-				edm4hep::ConstTrack oth = ( nhit0 > nhit1 ? trk1 :  trk0 ) ;
+				edm4hep::Track trk = ( nhit0 > nhit1 ? trk0 :  trk1 ) ;
+				edm4hep::Track oth = ( nhit0 > nhit1 ? trk1 :  trk0 ) ;
 
 				bool  outward = ( nhit0 > nhit1  ?  lthl0 <= lthf1 + overlapRows :  lthl1 <= lthf0 + overlapRows ) ;
 
@@ -529,7 +529,7 @@ namespace clupatra_new{
 			TrackCircleDistance(float dCut) : _dCutSquared( dCut*dCut ) , _dCut(dCut){}
 
 			/** Merge condition: ... */
-			bool operator() ( nnclu::Element<ClupaPlcioConstTrack>* h0, nnclu::Element<ClupaPlcioConstTrack>* h1);
+			bool operator() ( nnclu::Element<ClupaPlcioTrack>* h0, nnclu::Element<ClupaPlcioTrack>* h1);
   protected:
     float _dCutSquared ;
     float _dCut ;
@@ -537,7 +537,7 @@ namespace clupatra_new{
 	};
 
 	struct TrackZSort {  // sort tracks wtr to abs(z_average )
-		bool operator()( edm4hep::ConstTrack l, edm4hep::ConstTrack r);
+		bool operator()( edm4hep::Track l, edm4hep::Track r);
 	};
 
 
-- 
GitLab