diff --git a/Reconstruction/Tracking/CMakeLists.txt b/Reconstruction/Tracking/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3335b022e55e3dc1b442575d37c998c798ec33e9
--- /dev/null
+++ b/Reconstruction/Tracking/CMakeLists.txt
@@ -0,0 +1,27 @@
+gaudi_subdir(Tracking v0r0)
+
+find_package(GEAR REQUIRED)
+find_package(GSL REQUIRED )
+find_package(LCIO REQUIRED )
+find_package(EDM4HEP REQUIRED )
+find_package(DD4hep COMPONENTS DDCore DDRec REQUIRED)
+
+
+gaudi_depends_on_subdirs(
+ Service/GearSvc
+ Service/EventSeeder
+ Service/TrackSystemSvc
+)
+
+set(Tracking_srcs
+ src/Clupatra/*.cpp
+)
+
+# Modules
+gaudi_add_module(Tracking ${Tracking_srcs}
+ INCLUDE_DIRS GaudiKernel gear ${GSL_INCLUDE_DIRS} ${LCIO_INCLUDE_DIRS}
+ LINK_LIBRARIES GaudiAlgLib GaudiKernel ${GEAR_LIBRARIES} ${GSL_LIBRARIES} ${LCIO_LIBRARIES} TrackSystemSvcLib
+ -Wl,--no-as-needed
+ EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+ -Wl,--as-needed
+)
diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c048c09816dc2e85f628db1294cefbbb845f3bee
--- /dev/null
+++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
@@ -0,0 +1,1377 @@
+#include "ClupatraAlg.h"
+
+#include "clupatra_new.h"
+
+#include <time.h>
+#include <vector>
+#include <map>
+#include <set>
+#include <algorithm>
+#include <math.h>
+#include <cmath>
+#include <memory>
+#include <float.h>
+
+//---- MarlinUtil
+
+//---- LCIO ---
+#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"
+
+//-------gsl -----
+#include "gsl/gsl_randist.h"
+#include "gsl/gsl_cdf.h"
+
+
+//---- GEAR ----
+#include "gear/GEAR.h"
+#include "gear/TPCParameters.h"
+#include "gear/ZPlanarParameters.h"
+#include "gear/ZPlanarLayerLayout.h"
+#include "gear/PadRowLayout2D.h"
+#include "gear/BField.h"
+
+
+#include "TrackSystemSvc/IMarlinTrack.h"
+#include "TrackSystemSvc/IMarlinTrkSystem.h"
+#include "TrackSystemSvc/MarlinTrkUtils.h"
+
+#include "RuntimeMap.h"
+
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/TrackCollection.h"
+// #include "edm4hep/TrackerHitPlane.h"
+
+
+using namespace edm4hep ;
+using namespace MarlinTrk ;
+
+using namespace clupatra_new ;
+
+/*
+ namespace edm4hep::TrackState {
+ const int AtIP = 0;
+ const int AtFirstHit = 1;
+ const int AtLastHit = 2;
+ const int AtCalorimeter = 3;
+ };
+ */
+
+RuntimeMap<edm4hep::ConstTrack, clupatra_new::TrackInfoStruct*> TrackInfo_of_edm4hepTrack;
+RuntimeMap<edm4hep::ConstTrack, MarlinTrk::IMarlinTrack*> MarTrk_of_edm4hepTrack;
+RuntimeMap<MarlinTrk::IMarlinTrack*, clupatra_new::CluTrack*> CluTrk_of_MarTrack;
+RuntimeMap<edm4hep::ConstTrackerHit, clupatra_new::Hit*> GHitof;
+RuntimeMap<clupatra_new::CluTrack*, MarlinTrk::IMarlinTrack*> MarTrkof;
+
+gear::GearMgr* gearMgr;
+#define WRITE_PICKED_DEBUG_TRACKS false
+
+//----------------------------------------------------------------
+
+template <class In, class Pred> In find_smallest(In first, In last, Pred p, double& d){
+
+ In res = last ;
+
+ double min = DBL_MAX ;
+
+ while( first != last ){
+
+ double val = p( *first) ;
+
+ if( val < min ){
+
+ res = first ;
+ min = val ;
+ }
+
+ ++first ;
+ }
+
+ d = min ;
+
+ // FIXME: Mingrui ignore the debug
+
+ return res ;
+}
+//----------------------------------------------------------------
+struct Distance3D2{
+ gear::Vector3D _pos ;
+ Distance3D2( const gear::Vector3D& pos) : _pos( pos ) {}
+ template <class T>
+ double operator()( const T* t) {
+ gear::Vector3D p( t->getPosition() ) ;
+ return ( p - _pos ).r2() ;
+
+ }
+};
+
+//----------------------------------------------------------------
+struct StripDistance2{
+ gear::Vector3D _pos ;
+ StripDistance2( const gear::Vector3D& pos) : _pos( pos ) {}
+
+ double operator()( const TrackerHit* t) {
+
+/*
+ gear::Vector3D p(t->getPosition()[0], t->getPosition()[1], t->getPosition()[2]) ;
+
+ const TrackerHitPlane* h = (const TrackerHitPlane*) t ;
+
+ gear::Vector3D v( 1. , h->getU()[1] , h->getU()[0] , gear::Vector3D::spherical ) ;
+
+ double d = ( p - _pos ).dot(v) ;
+
+ return d*d ;
+*/
+return 1;
+ }
+};
+
+//----------------------------------------------------------------
+struct MeanAbsZOfTrack{
+ double operator()( ConstTrack t){
+ double z = 0 ;
+ int hitCount = 0 ;
+ /*
+ const TrackerHitVec& hV = t->getTrackerHits() ;
+ for(unsigned i=0, N = hV.size() ; i<N ; ++i){
+ z += hV[i]->getPosition()[2] ;
+ ++hitCount ;
+ }
+ */
+ for (auto iter = t.trackerHits_begin(); iter != t.trackerHits_end(); iter++) {
+ z += iter->getPosition()[2];
+ ++hitCount;
+ }
+ return ( hitCount>0 ? std::abs( z ) / hitCount : 0 ) ;
+ }
+};
+
+
+DECLARE_COMPONENT( ClupatraAlg )
+
+ClupatraAlg::ClupatraAlg(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc), _trksystem(0), _gearTPC(0) {
+
+ // _description = "ClupatraProcessor : nearest neighbour clustering seeded pattern recognition" ;
+
+ // Input Collections
+ // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+ declareProperty("TPCHitCollection", _TPCHitCollectionHandle, "handler of the tpc hit input collections");
+
+ // Output Collections
+ // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+ declareProperty("ClupatraTracks", _ClupatraTrackCollectionHandle, "handler of the collection with final TPC tracks");
+ declareProperty("ClupatraTrackSegments", _ClupatraTrackSegmentCollectionHandle, "handler of the collection that has the individual track segments");
+
+ }
+
+void ClupatraAlg::printParameters() {
+
+ debug() << "Print parameters:" << endmsg;
+ debug() << _distCut << endmsg;
+ debug() << _cosAlphaCut << endmsg;
+ debug() << _nLoop << endmsg;
+ debug() << _minCluSize << endmsg;
+ debug() << _duplicatePadRowFraction << endmsg;
+ debug() << _dChi2Max << endmsg;
+ debug() << _chi2Cut << endmsg;
+ debug() << _maxStep << endmsg;
+ debug() << _pickUpSiHits << endmsg;
+ debug() << _createDebugCollections << endmsg;
+ debug() << _padRowRange << endmsg;
+ debug() << _nZBins << endmsg;
+ debug() << _minLayerFractionWithMultiplicity << endmsg;
+ debug() << _minLayerNumberWithMultiplicity << endmsg;
+ debug() << _trackStartsInnerDist << endmsg;
+ debug() << _trackEndsOuterCentralDist << endmsg;
+ debug() << _trackEndsOuterForwardDist << endmsg;
+ debug() << _trackIsCurlerOmega << endmsg;
+
+}
+
+
+StatusCode ClupatraAlg::initialize() {
+
+ // Usually a good idea to
+ 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->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");
+
+ return GaudiAlgorithm::initialize();
+
+}
+
+
+StatusCode ClupatraAlg::execute() {
+
+ TrackInfo_of_edm4hepTrack.init();
+ MarTrk_of_edm4hepTrack.init();
+ CluTrk_of_MarTrack.init();
+ MarTrkof.init();
+ GHitof.init();
+
+
+ debug() << "Clupatra Algorithm started" << endmsg;
+
+ // clock_t start = clock() ;
+ Timer timer ;
+ unsigned t_init = timer.registerTimer(" initialization " ) ;
+ unsigned t_seedtracks = timer.registerTimer(" extend seed tracks " ) ;
+ unsigned t_recluster = timer.registerTimer(" recluster leftovers " ) ;
+ unsigned t_split = timer.registerTimer(" split clusters " ) ;
+ unsigned t_finalfit = timer.registerTimer(" final refit " ) ;
+ unsigned t_merge = timer.registerTimer(" merge segments " ) ;
+ unsigned t_pickup = timer.registerTimer(" pick up Si hits " ) ;
+
+ timer.start() ;
+
+ // the clupa wrapper hits that hold pointers to PLCIO hits plus some additional parameters
+ // create them in a vector for convenient memeory mgmt
+ std::vector<ClupaHit> clupaHits ;
+
+ // on top of the clupahits we need the tiny wrappers for clustering - they are created on the heap
+ // and we put them in a vector of pointers that takes ownership (i.e. deletes them at the end)
+ HitVec nncluHits ;
+ nncluHits.setOwner( true ) ;
+
+
+ // this is the final list of cluster tracks
+ Clusterer::cluster_list cluList ;
+ cluList.setOwner() ;
+
+
+ PLCIOTrackConverter converter ;
+
+ // Gear can be put as global
+ auto _gear = service<IGearSvc>("GearSvc");
+ gearMgr = _gear->getGearMgr();
+
+ _gearTPC = &(gearMgr->getTPCParameters());
+ _bfield = gearMgr->getBField().at( gear::Vector3D(0.,0.0,0.) ).z() ;
+
+ // Support for more than one module
+ // The ternary operator is used to make the trick with the static variable which
+ // is supposed to be calculated only once, also for performance reason
+ static const unsigned int maxTPCLayers = (gearMgr->getDetectorName() == "LPTPC" ) ?
+ _gearTPC->getModule(0).getNRows() + _gearTPC->getModule(2).getNRows() + _gearTPC->getModule(5).getNRows() : // LCTPC
+ _gearTPC->getModule(0).getNRows(); // ILD
+
+ double driftLength = _gearTPC->getMaxDriftLength() ;
+ ZIndex zIndex( -driftLength , driftLength , _nZBins ) ;
+
+
+ const edm4hep::TrackerHitCollection* col = nullptr;
+ debug() << "col" << endmsg;
+
+ try{ col = _TPCHitCollectionHandle.get();
+
+ } catch(...) {
+ // FIXME: Mingrui fix the output message
+ // streamlog_out( WARNING ) << " input collection not in event : " << _colName << " - nothing to do !!! " << std::endl ;
+ }
+
+ //===============================================================================================
+ // create clupa and clustering hits for every edm4hep hit
+ //===============================================================================================
+ if (col == nullptr) {
+ return StatusCode::SUCCESS;
+ }
+
+
+ int nHit = col->size() ;
+
+ 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.
+ for(int i=0 ; i < nHit ; ++i ) {
+
+ //------
+ ConstTrackerHit th(col->at(i));
+ // debug() << i << " " << th->getCellID0() << " " << th->getCellID1() << endmsg;
+ if ( fabs(th.getPosition()[2]) > driftLength ) continue;
+
+ ClupaHit* ch = & clupaHits[i] ;
+
+ Hit* gh = new Hit( ch ) ;
+
+ nncluHits.push_back( gh ) ;
+
+ //-------
+ // FIXME: Here we should have a resolution
+ GHitof(th) = gh ; // assign the clupa hit to the LCIO hit for memory mgmt
+
+ ch->edm4hepHit = th ;
+
+ ch->pos = gear::Vector3D(th.getPosition()[0], th.getPosition()[1], th.getPosition()[2]) ;
+
+ // int padIndex = padLayout.getNearestPad( ch->pos.rho() , ch->pos.phi() ) ;
+ // ch->layer = padLayout.getRowNumber( padIndex ) ;
+ // 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 ;
+
+ ch->zIndex = zIndex( &th ) ;
+
+ //ch->phiIndex = ....
+
+ }
+ debug() << "Get hits for clustering" << endmsg;
+
+ //---------------------------------------------------------------------------------------------------------
+
+ std::sort( nncluHits.begin(), nncluHits.end() , ZSort() ) ;
+
+ //---------------------------------------------------------------------------------------------------------
+
+ HitListVector hitsInLayer( maxTPCLayers ) ;
+ addToHitListVector( nncluHits.begin(), nncluHits.end() , hitsInLayer ) ;
+
+ // debug() << " added " << nncluHits.size() << " tp hitsInLayer - > size " << hitsInLayer.size() << std::endl ;
+
+ //---------------------------------------------------------------------------------------------------------
+
+ //===============================================================================================
+ // create output collections ( some optional )
+ //===============================================================================================
+
+ // FIXME: Mingrui how to create the correct output collection?
+
+ // const bool writeSeedCluster = _createDebugCollections ;
+ // const bool writeCluTrackSegments = _createDebugCollections ;
+ // const bool writeLeftoverClusters = _createDebugCollections ;
+ // const bool writeQualityTracks = _createDebugCollections ;
+ // const bool writeDebugTracks = WRITE_PICKED_DEBUG_TRACKS ;
+
+ static const bool copyTrackSegments = false ;
+
+ /*
+ * FIXME Mingrui remove these collections for now
+ LCCollectionVec* seedCol = ( writeSeedCluster ? newTrkCol( "ClupatraSeedCluster" , evt ) : 0 ) ;
+ LCCollectionVec* cluCol = ( writeCluTrackSegments ? newTrkCol( "ClupatraInitialTrackSegments" , evt ) : 0 ) ;
+ LCCollectionVec* locCol = ( writeCluTrackSegments ? newTrkCol( "ClupatraLeftoverClusters" , evt ) : 0 ) ;
+
+ LCCollectionVec* incSegCol = ( writeCluTrackSegments ? newTrkCol( "ClupatraIncompleteSegments" , evt , true ) : 0 ) ;
+ LCCollectionVec* curSegCol = ( writeCluTrackSegments ? newTrkCol( "ClupatraCurlerSegments" , evt , true ) : 0 ) ;
+ LCCollectionVec* finSegCol = ( writeCluTrackSegments ? newTrkCol( "ClupatraFinalTrackSegments" , evt , true ) : 0 ) ;
+
+ //LCCollectionVec* goodCol = ( writeQualityTracks ? newTrkCol( "ClupatraGoodQualityTracks" , evt ,true ) : 0 ) ;
+ //LCCollectionVec* fairCol = ( writeQualityTracks ? newTrkCol( "ClupatraFairQualityTracks" , evt ,true ) : 0 ) ;
+ LCCollectionVec* poorCol = ( writeQualityTracks ? newTrkCol( "ClupatraPoorQualityTracks" , evt , true ) : 0 ) ;
+
+ LCCollectionVec* debugCol= ( writeDebugTracks ? newTrkCol( "ClupatraDebugTracks" , evt , false ) : 0 ) ;
+ if( WRITE_PICKED_DEBUG_TRACKS )
+ DebugTracks::setCol( debugCol , this ) ;
+
+ LCCollectionVec* outerCol = ( _createDebugCollections ? newTrkCol( "ClupatraOuterSegments" , evt ,true ) : 0 ) ;
+ LCCollectionVec* innerCol = ( _createDebugCollections ? newTrkCol( "ClupatraInnerSegments" , evt ,true ) : 0 ) ;
+ LCCollectionVec* middleCol = ( _createDebugCollections ? newTrkCol( "ClupatraMiddleSegments" , evt ,true ) : 0 ) ;
+ */
+
+
+
+ TrackCollection* outCol = _ClupatraTrackCollectionHandle.createAndPut();
+ TrackCollection* tsCol = _ClupatraTrackSegmentCollectionHandle.createAndPut();
+ std::vector<ClupaPlcioConstTrack*> outCol_tmp;
+ std::vector<ClupaPlcioConstTrack*> tsCol_tmp;
+
+ //---------------------------------------------------------------------------------------------------------
+
+ timer.time(t_init ) ;
+
+ //===============================================================================================
+ // first main step of clupatra:
+ // * cluster in pad row range - starting from the outside - to find clean cluster segments
+ // * extend the track segments with best matching hits, based on extrapolation to next layer(s)
+ // * add the hits and apply a Kalman filter step ( track segement is always best estimate )
+ // * repeat in backward direction ( after smoothing back, to get a reasonable track
+ // state for extrapolating backwards )
+ //===============================================================================================
+
+ Clusterer nncl ;
+
+ int outerRow = 0 ;
+
+ nnclu::PtrVector<MarlinTrk::IMarlinTrack> seedTrks ;
+ seedTrks.setOwner() ; // memory mgmt - will delete MarlinTrks at the end
+
+ IMarlinTrkFitter fitter( _trksystem ) ;
+
+
+ debug() << "===============================================================================================\n"
+ << " first step of Clupatra algorithm: find seeds with NN-clustering in " << _nLoop << " loops - max dist = " << _distCut <<" \n"
+ << "===============================================================================================\n" << 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
+ //
+ double dcut = _distCut / _nLoop ;
+ for(int nloop=1 ; nloop <= _nLoop ; ++nloop) {
+
+ HitDistance dist( nloop * dcut , _cosAlphaCut ) ;
+
+ outerRow = maxTPCLayers - 1 ;
+
+ while( outerRow >= _minCluSize ) {
+
+ HitVec hits ;
+ hits.reserve( nHit ) ;
+
+ // add all hits in pad row range to hits
+ for(int iRow = outerRow ; iRow > ( outerRow - _padRowRange) ; --iRow ) {
+
+ if( iRow > -1 ) {
+
+ // debug() << " copy " << hitsInLayer[ iRow ].size() << " hits for row " << iRow << endmsg ;
+
+ std::copy( hitsInLayer[ iRow ].begin() , hitsInLayer[ iRow ].end() , std::back_inserter( hits ) ) ;
+ }
+ }
+
+ //----- cluster in given pad row range -----------------------------
+ Clusterer::cluster_list sclu ;
+ sclu.setOwner() ;
+
+ // FIXME Mingrui
+ //streamlog_out( DEBUG2 ) << " call cluster_sorted with " << hits.size() << " hits " << std::endl ;
+
+ nncl.cluster_sorted( hits.begin(), hits.end() , std::back_inserter( sclu ), dist , _minCluSize ) ;
+
+ const static int merge_seeds = true ;
+
+ if( merge_seeds ) { //-----------------------------------------------------------------------
+
+ // sometimes we have split seed clusters as one link is just above the cut
+ // -> recluster in all hits of small clusters with 1.2 * cut
+ float _smallClusterPadRowFraction = 0.9 ;
+ float _cutIncrease = 1.2 ;
+ // fixme: could make parameters ....
+
+ HitVec seedhits ;
+ Clusterer::cluster_list smallclu ;
+ smallclu.setOwner() ;
+ split_list( sclu, std::back_inserter(smallclu), ClusterSize( int( _padRowRange * _smallClusterPadRowFraction) ) ) ;
+ for( Clusterer::cluster_list::iterator sci=smallclu.begin(), end= smallclu.end() ; sci!=end; ++sci ){
+ for( Clusterer::cluster_type::iterator ci=(*sci)->begin(), end= (*sci)->end() ; ci!=end;++ci ){
+ seedhits.push_back( *ci ) ;
+ }
+ }
+ // free hits from bad clusters
+ std::for_each( smallclu.begin(), smallclu.end(), std::mem_fun( &CluTrack::freeElements ) ) ;
+
+ HitDistance distLarge( nloop * dcut * _cutIncrease ) ;
+
+ nncl.cluster_sorted( seedhits.begin(), seedhits.end() , std::back_inserter( sclu ), distLarge , _minCluSize ) ;
+
+ } //------------------------------------------------------------------------------------------
+
+ // FIXME Mingrui
+ // debug() << " found " << sclu.size() << " clusters " << std::endl ;
+
+ // try to split up clusters according to multiplicity
+ int layerWithMultiplicity = _padRowRange - 2 ; // fixme: make parameter
+ split_multiplicity( sclu , layerWithMultiplicity , 10 ) ;
+
+
+ // remove clusters whith too many duplicate hits per pad row
+ Clusterer::cluster_list bclu ; // bad clusters
+ bclu.setOwner() ;
+ split_list( sclu, std::back_inserter(bclu), DuplicatePadRows( maxTPCLayers, _duplicatePadRowFraction ) ) ;
+ // free hits from bad clusters
+ std::for_each( bclu.begin(), bclu.end(), std::mem_fun( &CluTrack::freeElements ) ) ;
+
+
+ // ---- now we also need to remove the hits from good cluster seeds from the hitsInLayers:
+ for( Clusterer::cluster_list::iterator sci=sclu.begin(), end= sclu.end() ; sci!=end; ++sci ){
+ for( Clusterer::cluster_type::iterator ci=(*sci)->begin(), end= (*sci)->end() ; ci!=end;++ci ){
+
+ // this is not cheap ...
+ hitsInLayer[ (*ci)->first->layer ].remove( *ci ) ;
+ }
+ }
+
+ // now we have 'clean' seed clusters
+ // Write debug collection with seed clusters:
+ // convert the clusters into tracks and write the resulting tracks into the existing debug collection seedCol.
+ // The conversion is performed by the STL transform() function, the insertion to the end of the
+ // debug track collection is done by creating an STL back_inserter iterator on the LCCollectionVector seedCol
+ /* Debug FIXME Mingrui
+ if( writeSeedCluster ) {
+ std::transform( sclu.begin(), sclu.end(), std::back_inserter( *seedCol ) , converter ) ;
+ }
+ */
+
+ // std::transform( sclu.begin(), sclu.end(), std::back_inserter( seedTrks) , fitter ) ;
+ // reduce memory footprint: deal with one KalTest track at a time and delete it, when done
+
+ for( Clusterer::cluster_list::iterator icv = sclu.begin(), end =sclu.end() ; icv != end ; ++ icv ) {
+ // icv = *cluster_list type
+
+ int nHitsAdded = 0 ;
+
+ // debug() << " call fitter for seed cluster with " << (*icv)->size() << " hits " << endmsg;
+ // int counter = 0;
+ // for( Clusterer::cluster_type::iterator ci=(*icv)->begin(), end= (*icv)->end() ; ci!=end; ++ci ) {
+ // debug() << counter++ << " " << *((*ci)->first->edm4hepHit) << " \nlayer " << (*ci)->first->layer << endmsg;
+ // }
+
+
+ MarlinTrk::IMarlinTrack* mTrk = fitter( *icv ) ;
+
+ // std::vector<std::pair<edm4hep::ConstTrackerHit, double> > hitsInFit ;
+ // mTrk->getHitsInFit( hitsInFit ) ;
+ // for (auto hit : hitsInFit) std::cout << hit.first << std::endl;
+
+ nHitsAdded += addHitsAndFilter( *icv , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex ) ;
+
+ static const bool backward = true ;
+ nHitsAdded += addHitsAndFilter( *icv , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex, backward ) ;
+ // in order to use smooth for backward extrapolation call with _trksystem - does not work well...
+ // nHitsAdded += addHitsAndFilter( *icv , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex, backward , _trksystem ) ;
+
+
+ // drop seed clusters with no hits added - but not in the very forward region...
+ // debug() << "Goes here" << endmsg;
+ if( nHitsAdded < 1 && outerRow > 2*_padRowRange ){ //FIXME: make parameter ?
+
+ ConstTrack edm4hepTrk( converter( *icv ) ) ;
+ // debug() << "Goes goes here" << endmsg;
+
+ // FIXME Mingrui
+ //streamlog_out( DEBUG2) << "============= poor seed cluster - no hits added - started from row " << outerRow << "\n"
+ //<< *edm4hepTrk << std::endl ;
+
+
+ for( Clusterer::cluster_type::iterator ci=(*icv)->begin(), end= (*icv)->end() ; ci!=end; ++ci ) {
+ hitsInLayer[ (*ci)->first->layer ].push_back( *ci ) ;
+ }
+ (*icv)->freeElements() ;
+ (*icv)->clear() ;
+ }
+
+ /*
+ if( writeCluTrackSegments ) // ---- store track segments from the first main step -----
+ cluCol->addElement( converter( *icv ) );
+ Debug FIXME Mingrui
+ */
+ // reset the pointer to the KalTest track - as we are done with this track
+ MarTrkof(*icv) = 0;
+
+ delete mTrk ;
+ }
+
+ // merge the good clusters to final list
+ cluList.merge( sclu ) ;
+
+ outerRow -= _padRowRange ;
+
+ } //while outerRow > padRowRange
+
+ } // nloop
+
+
+ timer.time( t_seedtracks ) ;
+
+ timer.time( t_recluster ) ;
+
+ //===============================================================================================
+ // do a global reclustering in leftover hits
+ //===============================================================================================
+ debug() << "do global reclustering" << endmsg;
+ static const int do_global_reclustering = true ;
+ if( do_global_reclustering ) {
+
+ outerRow = maxTPCLayers - 1 ;
+
+ int padRangeRecluster = 50 ; // FIXME: make parameter
+ // define an inner cylinder where we exclude hits from re-clustering:
+ double zMaxInnerHits = driftLength * .67 ; // FIXME: make parameter
+ double rhoMaxInnerHits = _gearTPC->getPlaneExtent()[0] + ( _gearTPC->getPlaneExtent()[1] - _gearTPC->getPlaneExtent()[0] ) * .67 ;// FIXME: make parameter
+
+
+ while( outerRow > 0 ) {
+
+
+ Clusterer::cluster_list loclu ; // leftover clusters
+ loclu.setOwner() ;
+
+ HitVec hits ;
+
+ int minRow = ( ( outerRow - padRangeRecluster ) > -1 ? ( outerRow - padRangeRecluster ) : -1 ) ;
+
+ // add all hits in pad row range to hits
+ for(int iRow = outerRow ; iRow > minRow ; --iRow ) {
+
+
+ for( HitList::iterator hlIt=hitsInLayer[ iRow ].begin() , end = hitsInLayer[ iRow ].end() ; hlIt != end ; ++hlIt ) {
+
+ if( std::abs( (*hlIt)->first->pos.z() ) > zMaxInnerHits || (*hlIt)->first->pos.rho() > rhoMaxInnerHits ) {
+ hits.push_back( *hlIt ) ;
+ }
+ }
+ }
+
+
+ HitDistance distSmall( _distCut ) ;
+ nncl.cluster( hits.begin(), hits.end() , std::back_inserter( loclu ), distSmall , _minCluSize ) ;
+
+ // Write debug collection using STL transform() function on the clusters
+ /* Debug FIXME Mingrui
+ if( writeLeftoverClusters )
+ std::transform( loclu.begin(), loclu.end(), std::back_inserter( *locCol ) , converter ) ;
+ */
+
+
+ // timer.time( t_recluster ) ;
+
+ //===============================================================================================
+ // now we split the clusters based on their hit multiplicities
+ //===============================================================================================
+
+
+ // _dChi2Max = 5. * _dChi2Max ; //FIXME !!!!!!!!!
+
+ for( Clusterer::cluster_list::iterator it= loclu.begin(), end= loclu.end() ; it != end ; ++it ){
+
+ CluTrack* clu = *it ;
+
+ // FIXME: Mingrui
+
+ std::vector<int> mult(8) ;
+ // get hit multiplicities up to 6 ( 7 means 7 or higher )
+ getHitMultiplicities( clu , mult ) ;
+
+ // FIXME Mingrui
+ // streamlog_out( DEBUG3 ) << " **** left over cluster with hit multiplicities: \n" ;
+ for( unsigned i=0,n=mult.size() ; i<n ; ++i) {
+ // FIXME Mingrui
+ // streamlog_out( DEBUG3 ) << " m["<<i<<"] = " << mult[i] << "\n" ;
+ }
+
+
+ if( float( mult[5]) / mult[0] >= _minLayerFractionWithMultiplicity && mult[5] > _minLayerNumberWithMultiplicity ) {
+
+ Clusterer::cluster_list reclu ; // reclustered leftover clusters
+ reclu.setOwner() ;
+
+ create_n_clusters( *clu , reclu , 5 ) ;
+
+ std::transform( reclu.begin(), reclu.end(), std::back_inserter( seedTrks) , fitter ) ;
+
+ for( Clusterer::cluster_list::iterator ir= reclu.begin(), end= reclu.end() ; ir != end ; ++ir ){
+
+ // FIXME Mingrui
+ // streamlog_out( DEBUG5 ) << " extending mult-5 clustre of length " << (*ir)->size() << std::endl ;
+
+ addHitsAndFilter( *ir , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex) ;
+ static const bool backward = true ;
+ addHitsAndFilter( *ir , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex, backward ) ;
+ }
+
+ cluList.merge( reclu ) ;
+ }
+
+ else if( float( mult[4]) / mult[0] >= _minLayerFractionWithMultiplicity && mult[4] > _minLayerNumberWithMultiplicity ) {
+
+ Clusterer::cluster_list reclu ; // reclustered leftover clusters
+ reclu.setOwner() ;
+
+ create_n_clusters( *clu , reclu , 4 ) ;
+
+ std::transform( reclu.begin(), reclu.end(), std::back_inserter( seedTrks) , fitter ) ;
+
+ for( Clusterer::cluster_list::iterator ir= reclu.begin(), end= reclu.end() ; ir != end ; ++ir ){
+
+
+ addHitsAndFilter( *ir , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex) ;
+ static const bool backward = true ;
+ addHitsAndFilter( *ir , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex, backward ) ;
+ }
+
+ cluList.merge( reclu ) ;
+ }
+
+ else if( float( mult[3]) / mult[0] >= _minLayerFractionWithMultiplicity && mult[3] > _minLayerNumberWithMultiplicity ) {
+
+ Clusterer::cluster_list reclu ; // reclustered leftover clusters
+ reclu.setOwner() ;
+
+ create_three_clusters( *clu , reclu ) ;
+
+ std::transform( reclu.begin(), reclu.end(), std::back_inserter( seedTrks) , fitter ) ;
+
+ for( Clusterer::cluster_list::iterator ir= reclu.begin(), end= reclu.end() ; ir != end ; ++ir ){
+
+ // FIXME Mingrui
+ // streamlog_out( DEBUG5 ) << " extending triplet clustre of length " << (*ir)->size() << std::endl ;
+
+ addHitsAndFilter( *ir , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex) ;
+ static const bool backward = true ;
+ addHitsAndFilter( *ir , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex, backward ) ;
+ }
+
+ cluList.merge( reclu ) ;
+ }
+
+ else if( float( mult[2]) / mult[0] >= _minLayerFractionWithMultiplicity && mult[2] > _minLayerNumberWithMultiplicity ) {
+
+ Clusterer::cluster_list reclu ; // reclustered leftover clusters
+ reclu.setOwner() ;
+
+ create_two_clusters( *clu , reclu ) ;
+
+ std::transform( reclu.begin(), reclu.end(), std::back_inserter( seedTrks) , fitter ) ;
+
+ for( Clusterer::cluster_list::iterator ir= reclu.begin(), end= reclu.end() ; ir != end ; ++ir ){
+
+ // FIXME Mingrui
+ // streamlog_out( DEBUG5 ) << " extending doublet clustre of length " << (*ir)->size() << std::endl ;
+
+ addHitsAndFilter( *ir , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex) ;
+ static const bool backward = true ;
+ addHitsAndFilter( *ir , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex, backward ) ;
+ }
+
+ cluList.merge( reclu ) ;
+
+ }
+ else if( float( mult[1]) / mult[0] >= _minLayerFractionWithMultiplicity && mult[1] > _minLayerNumberWithMultiplicity ) {
+
+
+ seedTrks.push_back( fitter( *it ) );
+
+ addHitsAndFilter( *it , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex) ;
+ static const bool backward = true ;
+ addHitsAndFilter( *it , hitsInLayer , _dChi2Max, _chi2Cut , _maxStep , zIndex, backward ) ;
+
+ cluList.push_back( *it ) ;
+
+ it = loclu.erase( it ) ;
+ --it ; // erase returns iterator to next element
+
+ } else {
+
+ // discard cluster and free hits
+ clu->freeElements() ;
+ }
+
+ }
+
+
+ outerRow -= padRangeRecluster ;
+
+ }
+ }
+
+ //=======================================================================================================================
+ // try again to gobble up hits at the ends .... - does not work right now, as there are no fits for the clusters....
+ //=======================================================================================================================
+
+
+ timer.time( t_split ) ;
+
+ //===============================================================================================
+ // now refit the tracks
+ //===============================================================================================
+
+ // FIXME Mingrui
+ debug() << " =========== refitting final " << cluList.size() << " track segments " << endmsg ;
+
+ //---- refit cluster tracks individually to save memory ( KalTest tracks have ~1MByte each)
+
+ IMarlinTrkFitter fit(_trksystem, _dChi2Max) ; // fixme: do we need a different chi2 max here ????
+
+ for( Clusterer::cluster_list::iterator icv = cluList.begin() , end = cluList.end() ; icv != end ; ++ icv ) {
+
+ if( (*icv)->empty() )
+ continue ;
+
+ MarlinTrk::IMarlinTrack* trk = fit( *icv ) ;
+ trk->smooth() ;
+ edm4hep::ConstTrack edm4hepTrk = converter( *icv ) ;
+ tsCol_tmp.push_back( new ClupaPlcioConstTrack(edm4hepTrk) ) ;
+ MarTrk_of_edm4hepTrack(edm4hepTrk) = 0 ;
+ delete trk ;
+ }
+
+ timer.time( t_finalfit) ;
+
+ //===============================================================================================
+ // optionally create collections of used and unused TPC hits
+ //===============================================================================================
+
+
+ //===============================================================================================
+ // compute some track parameters for possible merging
+ //===============================================================================================
+
+ typedef nnclu::NNClusterer<ClupaPlcioConstTrack> TrackClusterer ;
+ TrackClusterer nntrkclu ;
+ MakePLCIOElement<ClupaPlcioConstTrack> trkMakeElement ;
+
+ for( int i=0,N=tsCol_tmp.size() ; i<N ; ++i ) {
+
+ ConstTrack track = tsCol_tmp.at(i)->edm4hepTrack;
+ computeTrackInfo(track) ;
+ }
+
+ //===============================================================================================
+ // merge split segements
+ //===============================================================================================
+
+
+ static const int merge_split_segments = true ;
+
+ if( merge_split_segments ) {
+
+ for(unsigned l=0 ; l < 2 ; ++l ) { // do this twice ....
+
+ int nMax = tsCol_tmp.size() ;
+
+ TrackClusterer::element_vector incSegVec ;
+ incSegVec.setOwner() ;
+ incSegVec.reserve( nMax ) ;
+ TrackClusterer::cluster_vector incSegCluVec ;
+ incSegCluVec.setOwner() ;
+
+ for( int i=0,N=tsCol_tmp.size() ; i<N ; ++i ){
+
+ edm4hep::ConstTrack trk = tsCol_tmp.at(i)->edm4hepTrack;
+
+ const TrackInfoStruct* ti = TrackInfo_of_edm4hepTrack(trk);
+
+ bool isIncompleteSegment = !ti->isCurler && ( !ti->startsInner || ( !ti->isCentral && !ti->isForward ) ) ;
+
+ std::bitset<32> type = trk.getType() ;
+
+ if( isIncompleteSegment && ! type[ lcio::ILDTrackTypeBit::SEGMENT ]){
+
+ incSegVec.push_back( trkMakeElement( tsCol_tmp[i] ) ) ;
+
+ // FIXME Mingrui Debug part
+ /*
+ if( writeCluTrackSegments ) incSegCol->addElement( trk ) ;
+ */
+ }
+ }
+
+
+ TrackSegmentMerger trkMerge( _dChi2Max , _trksystem , _bfield ) ;
+
+ nntrkclu.cluster( incSegVec.begin() , incSegVec.end() , std::back_inserter( incSegCluVec ), trkMerge , 2 ) ;
+
+ // FIXME: Mingrui
+ // streamlog_out( DEBUG4 ) << " ===== merged track segments - # cluster: " << incSegCluVec.size()
+ // << " from " << incSegVec.size() << " incomplete track segments " << " ============================== " << std::endl ;
+
+ for( TrackClusterer::cluster_vector::iterator it= incSegCluVec.begin() ; it != incSegCluVec.end() ; ++it) {
+
+ // FIXME: Mingrui
+ // streamlog_out( DEBUG4 ) << edm4hep::header<edm4hep::Track>() << std::endl ;
+
+ TrackClusterer::cluster_type* incSegClu = *it ;
+
+ std::vector<edm4hep::ConstTrack> mergedTrk ;
+
+ // vector to collect hits from segments
+ // std::vector< TrackerHit* > hits ;
+ // hits.reserve( 1024 ) ;
+ // IMPL::TrackImpl* track = new IMPL::TrackImpl ;
+ // tsCol->addElement( track ) ;
+
+ CluTrack hits ;
+
+ for( TrackClusterer::cluster_type::iterator itC = incSegClu->begin() ; itC != incSegClu->end() ; ++ itC ){
+
+ //streamlog_out( DEBUG3 ) << lcshort( (*itC)->first ) << std::endl ;
+
+ edm4hep::ConstTrack trk = (*itC)->first->edm4hepTrack ;
+
+ mergedTrk.push_back( trk ) ;
+
+ // std::copy( trk->getTrackerHits().begin() , trk->getTrackerHits().end() , std::back_inserter( hits ) ) ;
+
+ /*
+ for( edm4hep::TrackerHitVec::const_iterator it = trk->getTrackerHits().begin() , END = trk->getTrackerHits().end() ; it != END ; ++ it ){
+ hits.addElement( GHitof(*it) ) ;
+ }
+ */
+ for (auto it = trk.trackerHits_begin(); it != trk.trackerHits_end(); it++) {
+ ConstTrackerHit hit = *it;
+ hits.addElement( GHitof(hit) );
+ }
+
+ // flag the segments so they can be ignored for final list
+ // FIXME I need to set type
+ // trk->setType( lcio::ILDTrackTypeBit::SEGMENT ) ;
+
+ // add old segments to new track
+ // track->addTrack( trk ) ;
+ }
+
+ // MarlinTrk::IMarlinTrack* mTrk = _trksystem->createTrack();
+ // EVENT::FloatVec icov( 15 ) ;
+ // icov[ 0] = 1e2 ;
+ // icov[ 2] = 1e2 ;
+ // icov[ 5] = 1e2 ;
+ // icov[ 9] = 1e2 ;
+ // icov[14] = 1e2 ;
+ // int result = createFinalisedLCIOTrack( mTrk, hits, track, !MarlinTrk::IMarlinTrack::backward, icov, _bfield, _dChi2Max ) ;
+ // //int result = createFinalisedLCIOTrack( mTrk, hits, track, ! MarlinTrk::IMarlinTrack::backward, icov, _bfield, _dChi2Max ) ;
+ // // ???
+
+ MarlinTrk::IMarlinTrack* mTrk = fit( &hits ) ;
+ mTrk->smooth() ;
+ edm4hep::ConstTrack track = converter( &hits ) ;
+ tsCol_tmp.push_back( new ClupaPlcioConstTrack(track) ) ;
+ MarTrk_of_edm4hepTrack(track) = 0 ;
+ delete mTrk ;
+ computeTrackInfo( track ) ;
+
+ // FIXME: Mingrui
+ // streamlog_out( DEBUG4 ) << " ****** created new track : " << " : " << lcshort( (Track*) track ) << std::endl ;
+
+ }
+ }// loop over l
+ }
+ //===============================================================================================
+ // merge curler segments
+ //===============================================================================================
+
+
+ static const int merge_curler_segments = true ;
+
+ if( merge_curler_segments ) {
+
+ int nMax = tsCol_tmp.size() ;
+
+ TrackClusterer::element_vector curSegVec ;
+ curSegVec.setOwner() ;
+ curSegVec.reserve( nMax ) ;
+ TrackClusterer::cluster_vector curSegCluVec ;
+ curSegCluVec.setOwner() ;
+
+
+ for( int i=tsCol_tmp.size()-1 ; i>=0 ; --i ){
+
+ ConstTrack trk = tsCol_tmp.at(i)->edm4hepTrack;
+
+ std::bitset<32> type = trk.getType() ;
+
+ if( type[ lcio::ILDTrackTypeBit::SEGMENT ] )
+ continue ; // ignore previously merged track segments
+
+ const TrackInfoStruct* ti = TrackInfo_of_edm4hepTrack(trk) ;
+
+ bool isCompleteTrack = ti && !ti->isCurler && ( ti->startsInner && ( ti->isCentral || ti->isForward ) );
+
+ if( !isCompleteTrack ){
+
+ curSegVec.push_back( trkMakeElement( tsCol_tmp[i] ) ) ;
+
+ // FIXME
+ // if( writeCluTrackSegments ) curSegCol->push_back( trk ) ;
+
+ } else { // ... is not a curler -> add a copy to the final tracks collection
+
+
+ if( copyTrackSegments) {
+
+ outCol_tmp.push_back( new ClupaPlcioConstTrack(trk) ) ;
+
+ }else{
+
+ outCol_tmp.push_back( new ClupaPlcioConstTrack(trk) ) ;
+
+ // FIXME !!! Possibly different with the LCIO version
+ // tsCol-> Do not remove it at the Beginning
+ // Very very important
+ // tsCol->removeElementAt( i ) ;
+ }
+
+ // FIXME Mingrui remove
+ // if( writeCluTrackSegments ) finSegCol->push_back( trk ) ;
+ }
+ }
+
+ //======================================================================================================
+
+
+ TrackCircleDistance trkMerge( 0.1 ) ;
+
+ 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;
+
+// debug() << "Finish this line" << endmsg;
+
+ for( TrackClusterer::cluster_vector::iterator it= curSegCluVec.begin() ; it != curSegCluVec.end() ; ++it) {
+
+ // FIXME: Mingrui
+ // streamlog_out( DEBUG4 ) << edm4hep::header<Track>() << std::endl ;
+
+ TrackClusterer::cluster_type* curSegClu = *it ;
+
+ std::list<ConstTrack> 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 ) ;
+ }
+
+
+ mergedTrk.sort( TrackZSort() ) ;
+
+ //================================================================================
+
+
+ if( copyTrackSegments) {
+
+ // ====== create a new LCIO track for the merged cluster ...
+ // edm4hep::Track trk;
+ edm4hep::Track trk;
+
+ trk.setType( lcio::ILDDetID::TPC ) ;
+
+ // == and copy all the hits
+ unsigned hitCount = 0 ;
+ for( std::list<ConstTrack>::iterator itML = mergedTrk.begin() ; itML != mergedTrk.end() ; ++ itML ){
+
+ for (auto itHit = (*itML).trackerHits_begin(); itHit != (*itML).trackerHits_end(); itHit++) {
+ trk.addToTrackerHits(*itHit);
+ }
+ hitCount += (*itML).trackerHits_size() ;
+
+ // add a pointer to the original track segment
+ trk.addToTracks(*itML) ;
+ }
+
+ // take track states from first and last track :
+ ConstTrack firstTrk = mergedTrk.front() ;
+ ConstTrack lastTrk = mergedTrk.back() ;
+
+ /* !!!!!!!!!!!!!! critical important FIXME should wait zoujiaheng
+ edm4hep::TrackState ts;
+ ts = firstTrk->getTrackState( lcio::TrackState::AtIP ) ;
+ if( ts ) trk->addTrackState( new TrackState( *ts ) ) ;
+
+ ts = firstTrk->getTrackState( lcio::TrackState::AtFirstHit ) ;
+ if( ts ) trk->addTrackState( new TrackState( *ts ) ) ;
+
+ ts = lastTrk->getTrackState( lcio::TrackState::AtLastHit ) ;
+ if( ts ) trk->addTrackState( new TrackState( *ts ) ) ;
+
+ ts = lastTrk->getTrackState( lcio::TrackState::AtCalorimeter ) ;
+ if( ts ) trk->addTrackState( new TrackState( *ts ) ) ;
+ */
+ if (hasTrackStateAt(firstTrk, lcio::TrackState::AtIP )) trk.addToTrackStates(getTrackStateAt(firstTrk, lcio::TrackState::AtIP));
+ if (hasTrackStateAt(firstTrk, lcio::TrackState::AtFirstHit )) trk.addToTrackStates(getTrackStateAt(firstTrk, lcio::TrackState::AtFirstHit));
+ if (hasTrackStateAt(lastTrk, lcio::TrackState::AtLastHit )) trk.addToTrackStates(getTrackStateAt(lastTrk, lcio::TrackState::AtLastHit));
+ if (hasTrackStateAt(lastTrk, lcio::TrackState::AtCalorimeter )) trk.addToTrackStates(getTrackStateAt(lastTrk, lcio::TrackState::AtCalorimeter));
+
+
+ MarTrk_of_edm4hepTrack(trk) = MarTrk_of_edm4hepTrack(firstTrk) ;
+
+ // FIXME: Mingrui Maybe this info not needed
+ //int hitsInFit = firstTrk->getSubdetectorHitNumbers()[ 2 * lcio::ILDDetID::TPC - 1 ] ;
+ trk.setChi2( firstTrk.getChi2() ) ;
+ trk.setNdf( firstTrk.getNdf() ) ;
+ trk.setDEdx( firstTrk.getDEdx() ) ;
+ trk.setDEdxError(firstTrk.getDEdxError()) ;
+
+ /* FIXME: Mingrui Maybe this subdetectorHitNumbers is not needed?
+ trk->subdetectorHitNumbers().resize( 2 * lcio::ILDDetID::ETD ) ;
+ trk->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::TPC - 2 ] = hitsInFit ;
+ trk->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::TPC - 1 ] = hitCount ;
+ */
+
+ double RMin = 0.0;
+ if (hasTrackStateAt(trk, lcio::TrackState::AtFirstHit ) ) {
+ TrackState ts = getTrackStateAt(trk, lcio::TrackState::AtFirstHit ) ;
+ RMin = sqrt( ts.referencePoint[0] * ts.referencePoint[0]
+ + ts.referencePoint[1] * ts.referencePoint[1] );
+
+ }
+ /*
+ double RMin = ( ts ?
+ sqrt( ts->referencePoint[0] * ts->referencePoint[0]
+ + ts->referencePoint[1] * ts->referencePoint[1] )
+ : 0.0 ) ;
+ */
+ trk.setRadiusOfInnermostHit( RMin ) ;
+
+
+
+ // 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 ;
+
+
+ outCol_tmp.push_back( new ClupaPlcioConstTrack(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() ;
+
+ ConstTrack trk = (ConstTrack) *itML++ ;
+
+ for( ; itML != mergedTrk.end() ; ++itML ){
+
+ // add a pointer to the original track segment
+ // FIXME I need to add sub track
+ // trk->addTrack(**itML) ;
+ }
+
+ outCol_tmp.push_back( new ClupaPlcioConstTrack(trk) ) ;
+
+ //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 ;
+ }
+ }
+
+ }//================================================================================
+
+ }
+ //---------------------------------------------------------------------------------------------
+ // // add all tracks that have not been merged :
+
+// 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;
+ if( (*it)->second == 0 ){
+
+ // std::cout << "Can I get the track? 2" << std::endl;
+ // std::cout << (*it)->first->edm4hepTrack;
+ ConstTrack trk = (*it)->first->edm4hepTrack ;
+ // std::cout << "Can I get the track? 3" << std::endl;
+ // std::cout << trk << std::endl;
+
+ if( copyTrackSegments) {
+
+ ConstTrack t = trk;
+
+ TrackInfo_of_edm4hepTrack(t) = 0 ; // set extension to 0 to prevent double free ...
+
+ MarTrk_of_edm4hepTrack(t) = 0 ; // dynamic_cast<Track*>( (*it)->first )->ext<MarTrk>() ;
+
+ // 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) ) ;
+
+ } else {
+ outCol_tmp.push_back( new ClupaPlcioConstTrack(trk) ) ;
+
+ //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->removeElementAt( i ) ;
+ break ;
+ }
+ }
+ }
+
+
+ }
+ }
+ }
+
+ timer.time( t_merge ) ;
+
+// debug() << "Finish merging works here" << endmsg;
+
+ //---------------------------------------------------------------------------------------------------------
+ // pick up hits from Si trackers
+ //---------------------------------------------------------------------------------------------------------
+
+ /*
+ if( _pickUpSiHits ){
+ pickUpSiTrackerHits(outCol) ;
+ }
+ */
+ //---------------------------------------------------------------------------------------------------------
+
+ timer.time( t_pickup ) ;
+
+ totalCandidates = outCol_tmp.size();
+ 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);
+ }
+ for (auto t : outCol_tmp) {
+ delete t;
+ }
+ for (auto t : tsCol_tmp) {
+ delete t;
+ }
+
+ _nEvt++ ;
+
+ return StatusCode::SUCCESS;
+}
+
+// ####### 001
+/*************************************************************************************************/
+
+void ClupatraAlg::computeTrackInfo( edm4hep::ConstTrack lTrk ){
+
+ if( ! TrackInfo_of_edm4hepTrack(lTrk) )
+ TrackInfo_of_edm4hepTrack(lTrk) = new TrackInfoStruct ;
+
+ float r_inner = _gearTPC->getPlaneExtent()[0] ;
+ float r_outer = _gearTPC->getPlaneExtent()[1] ;
+ float driftLength = _gearTPC->getMaxDriftLength() ;
+
+ // compute z-extend of this track segment
+ // auto hv = lTrk->getTrackerHits() ;
+
+ float zMin = 1e99 ;
+ float zMax = -1e99 ;
+ float zAvg = 0. ;
+
+ if( lTrk.trackerHits_size() > 1 ) {
+ zMin = lTrk.getTrackerHits( 0 ).getPosition()[2] ;
+ zMax = lTrk.getTrackerHits( lTrk.trackerHits_size() -1 ).getPosition()[2] ;
+ zAvg = ( zMax + zMin ) / 2. ;
+ }
+
+ if( zMin > zMax ){ // swap
+ float d = zMax ;
+ zMax = zMin ;
+ zMin = d ;
+ }
+
+
+ // protect against bad tracks
+ // if( tsF == 0 ) return ;
+ // if( tsL == 0 ) return ;
+ if (! hasTrackStateAt(lTrk, lcio::TrackState::AtFirstHit)) return; // StatusCode::FAILURE;
+ if (! hasTrackStateAt(lTrk, lcio::TrackState::AtFirstHit)) return; // StatueCode::FAILURE;
+
+ edm4hep::TrackState tsF = getTrackStateAt(lTrk, lcio::TrackState::AtFirstHit ) ;
+ edm4hep::TrackState tsL = getTrackStateAt(lTrk, lcio::TrackState::AtLastHit ) ;
+
+ gear::Vector3D fhPos(tsF.referencePoint[0], tsF.referencePoint[1], tsF.referencePoint[2]) ;
+ gear::Vector3D lhPos(tsL.referencePoint[0], tsL.referencePoint[1], tsL.referencePoint[2]) ;
+
+
+ TrackInfoStruct* ti = TrackInfo_of_edm4hepTrack(lTrk);
+
+ ti->startsInner = std::abs( fhPos.rho() - r_inner ) < _trackStartsInnerDist ; // first hit close to inner field cage
+ ti->isCentral = std::abs( lhPos.rho() - r_outer ) < _trackEndsOuterCentralDist ; // last hit close to outer field cage
+ ti->isForward = driftLength - std::abs( lhPos.z() ) < _trackEndsOuterForwardDist ; // last hitclose to endcap
+ ti->isCurler = std::abs( tsF.omega ) > _trackIsCurlerOmega ; // curler segment ( r <~ 1m )
+
+ ti->zMin = zMin ;
+ ti->zMax = zMax ;
+ ti->zAvg = zAvg ;
+
+}
+
+
+/*************************************************************************************************/
+// void ClupatraAlg::check() {
+
+
+// check that all Clupatra tracks actually have the four canonical track states set:
+
+// FIXME debug
+// streamlog_out( DEBUG5 ) << " ------------------- check() called " << std::endl ;
+
+// for( LCIterator<Track> it( evt, _segmentsOutColName ) ; EVENT::Track* trk = it.next() ; ) {
+/* FIXME debug
+ for( LCIterator<Track> it( evt, _outColName ) ; EVENT::Track* trk = it.next() ; ) {
+
+ const EVENT::TrackState* ts0 = trk->getTrackState( edm4hep::TrackState::AtIP ) ;
+ const EVENT::TrackState* ts1 = trk->getTrackState( edm4hep::TrackState::AtFirstHit ) ;
+ const EVENT::TrackState* ts2 = trk->getTrackState( edm4hep::TrackState::AtLastHit ) ;
+ const EVENT::TrackState* ts3 = trk->getTrackState( edm4hep::TrackState::AtCalorimeter ) ;
+
+
+// streamlog_out( DEBUG2 ) << lcshort( trk ) << ", " << ts0 << ", " << ts1 << ", " << ts2 << ", " << ts3 << std::endl ;
+// FIXME debug
+// streamlog_out( DEBUG3 ) << " -- " << ts0 << ", " << ts1 << ", " << ts2 << ", " << ts3 << std::endl ;
+
+if( ! ts0 || ! ts1 || ! ts2 || ! ts3 )
+
+// FIXME debug
+// streamlog_out( ERROR ) << " clupatra track w/ missing track state : " << lcshort( trk )
+// << " ts0-ts3 : " << ts0 << ", " << ts1 << ", " << ts2 << ", " << ts3 << std::endl ;
+
+
+}
+
+streamlog_out( DEBUG5 ) << " ------------------- check() done " << std::endl ;
+*/
+
+/*************************************************************************************************/
+// }
+
+
+
+//====================================================================================================
+
+StatusCode ClupatraAlg::finalize(){
+ tree->SaveAs("Tuple.root");
+
+ /*
+ * FIXME debug
+ streamlog_out( MESSAGE ) << "ClupatraProcessor::end() " << name()
+ << " processed " << _nEvt << " events in " << _nRun << " runs "
+ << std::endl ;
+ */
+ return StatusCode::SUCCESS;
+
+}
+
+
+//====================================================================================================
diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..80c10b64fb45156095a2ab13eec8a26783081266
--- /dev/null
+++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
@@ -0,0 +1,174 @@
+#ifndef ClupatraProcessor_h
+#define ClupatraProcessor_h 1
+
+#include "assert.h"
+#include "FWCore/DataHandle.h"
+#include "GearSvc/IGearSvc.h"
+#include "TrackSystemSvc/ITrackSystemSvc.h"
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/Track.h"
+#include "edm4hep/TrackerHitConst.h"
+#include "edm4hep/TrackCollection.h"
+#include <string>
+
+#include "gear/TPCModule.h"
+#include "RuntimeMap.h"
+#include "clupatra_new.h"
+#include "TrackingHelper.h"
+
+#include "TrackSystemSvc/MarlinTrkUtils.h"
+#include "TrackSystemSvc/HelixTrack.h"
+#include "TrackSystemSvc/HelixFit.h"
+#include "TrackSystemSvc/IMarlinTrack.h"
+
+namespace MarlinTrk{
+ class IMarlinTrkSystem ;
+}
+
+namespace gear{
+ class TPCParameters ;
+ class PadRowLayout2D ;
+}
+
+/** ClupatraAlgorithm : nearest neighbour clustering seeded pattern recognition using Kalman-Filter for extrapolation and
+ * hit search.
+ *
+ * @parameter TPCHitCollection Name of the tpc hit input collections
+ * @parameter OutputCollection Name of the output collection with final TPC tracks
+ * @parameter SegmentCollectionName Name of the output collection that has the individual track segments
+ * @parameter CreateDebugCollections optionally create some debug collection with intermediate track segments and used and unused hits
+ *
+ * @parameter DistanceCut Cut for distance between hits in mm for the seed finding
+ * @parameter Chi2Cut the maximum chi2-distance for which a hit is considered for merging
+ * @parameter CosAlphaCut Cut for max.angle between hits in consecutive layers for seed finding - NB value should be smaller than 1 - default is 0.9999999
+ * @parameter MaxDeltaChi2 the maximum delta Chi2 after filtering for which a hit is added to a track segement
+ *
+ * @parameter DuplicatePadRowFraction allowed fraction of hits in same pad row per track
+ * @parameter NLoopForSeeding number of seed finding loops - every loop increases the distance cut by DistanceCut/NLoopForSeeding
+ * @parameter NumberOfZBins number of bins in z over total length of TPC - hits from different z bins are nver merged
+ * @parameter PadRowRange number of pad rows used in initial seed clustering
+ *
+ * @parameter MaxStepWithoutHit the maximum number of layers without finding a hit before hit search search is stopped
+ * @parameter MinLayerFractionWithMultiplicity minimum fraction of layers that have a given multiplicity, when forcing a cluster into sub clusters
+ * @parameter MinLayerNumberWithMultiplicity minimum number of layers that have a given multiplicity, when forcing a cluster into sub clusters
+ * @parameter MinimumClusterSize minimum number of hits per cluster
+ *
+ * @parameter TrackEndsOuterCentralDist maximum radial distance [mm] from outer field cage of last hit, such that the track is considered to end at the end
+ * @parameter TrackEndsOuterForwardDist maximum distance in z [mm] from endplate of last hit, such that the track is considered to end at the end
+ * @parameter TrackIsCurlerOmega minimum curvature omega of a track segment for being considered a curler
+ * @parameter TrackStartsInnerDist maximum radial distance [mm] from inner field cage of first hit, such that the track is considered to start at the beginning
+ *
+ * @parameter EnergyLossOn Use Energy Loss in Fit
+ * @parameter MultipleScatteringOn Use MultipleScattering in Fit
+ * @parameter SmoothOn Smooth All Mesurement Sites in Fit
+ *
+ * @parameter pickUpSiHits try to pick up hits from Si-trackers
+ * @parameter SITHitCollection name of the SIT hit collections - used to extend TPC tracks if (pickUpSiHits==true)
+ * @parameter VXDHitCollection name of the VXD hit collections - used to extend TPC tracks if (pickUpSiHits==true)
+ *
+ * @parameter Verbosity verbosity level of this processor ("DEBUG0-4,MESSAGE0-4,WARNING0-4,ERROR0-4,SILENT")
+ *
+ * @author F.Gaede, DESY, 2011/2012
+ * @version $Id: ClupatraProcessor.h 4488 2013-04-05 12:03:59Z volynets $
+ */
+class ClupatraAlg : public GaudiAlgorithm {
+ // friend class AlgFactory<ClupatraAlg>;
+ public:
+
+
+ ClupatraAlg(const std::string& name, ISvcLocator* svcLoc) ;
+
+ /** Called at the begin of the job before anything is read.
+ * Use to initialize the processor, e.g. book histograms.
+ */
+ virtual StatusCode initialize() ;
+
+ /** Called for every run.
+ */
+ // virtual void processRunHeader( LCRunHeader* run ) ;
+
+ /** Called for every event - the working horse.
+ */
+ virtual StatusCode execute() ;
+
+ // virtual void check( LCEvent * evt ) ;
+
+ /** Called after data processing for clean up.
+ */
+ virtual StatusCode finalize() ;
+
+
+ protected:
+
+ /** helper method to compute a few track segment parameters (start and end points, z spread,...)
+ */
+ void computeTrackInfo(edm4hep::ConstTrack lTrk) ;
+
+
+ StatusCode pickUpSiTrackerHits(edm4hep::TrackCollection* trackCol) ;
+ void printParameters();
+
+ /** Input collection name.
+ */
+ std::string _colName ;
+ std::string _vxdColName ;
+ std::string _sitColName ;
+ std::string _outColName ;
+ std::string _segmentsOutColName ;
+
+ Gaudi::Property<float> _distCut{this, "DistanceCut", 40.0}; // Cut for distance between hits in mm for the seed finding
+ Gaudi::Property<float> _cosAlphaCut{this, "CosAlphaCut", 0.9999999}; // Cut for max.angle between hits in consecutive layers for seed finding - NB value should be smaller than 1 - default is 0.9999999 !!!
+ Gaudi::Property<int> _nLoop{this, "NLoopForSeeding", 4}; // number of seed finding loops - every loop increases the distance cut by DistanceCut/NLoopForSeeding
+ Gaudi::Property<int> _minCluSize{this, "MinimumClusterSize", 6}; // minimum number of hits per cluster
+ Gaudi::Property<float> _duplicatePadRowFraction{this, "DuplicatePadRowFraction", 0.1}; // allowed fraction of hits in same pad row per track
+ Gaudi::Property<float> _dChi2Max{this, "MaxDeltaChi2", 35.}; // the maximum delta Chi2 after filtering for which a hit is added to a track segement
+ Gaudi::Property<float> _chi2Cut{this, "Chi2Cut", 100.}; // the maximum chi2-distance for which a hit is considered for merging
+ Gaudi::Property<int> _maxStep{this, "MaxStepWithoutHit", 6}; // the maximum number of layers without finding a hit before hit search search is stopped
+ Gaudi::Property<bool> _pickUpSiHits{this, "pickUpSiHits", false}; // "try to pick up hits from Si-trackers"
+ Gaudi::Property<bool> _createDebugCollections{this, "CreateDebugCollections", false}; // optionally create some debug collection with intermediate track segments and used and unused hits
+ Gaudi::Property<int> _padRowRange{this, "PadRowRange", 15}; // "number of pad rows used in initial seed clustering"
+ Gaudi::Property<int> _nZBins{this, "NumberOfZBins", 150}; // "number of bins in z over total length of TPC - hits from different z bins are nver merged"
+ Gaudi::Property<float> _minLayerFractionWithMultiplicity{this, "MinLayerFractionWithMultiplicity" , 0.5}; // "minimum fraction of layers that have a given multiplicity, when forcing a cluster into sub clusters"
+ Gaudi::Property<float> _minLayerNumberWithMultiplicity{this, "MinLayerNumberWithMultiplicity", 3}; // "minimum number of layers that have a given multiplicity, when forcing a cluster into sub clusters"
+ Gaudi::Property<float> _trackStartsInnerDist{this, "TrackStartsInnerDist", (float) 25.}; // "maximum radial distance [mm] from inner field cage of first hit, such that the track is considered to start at the beginning "
+ 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"
+
+
+
+ DataHandle<edm4hep::TrackerHitCollection> _TPCHitCollectionHandle{"TPCTrackerHits", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> _SITHitCollectionHandle{"SIDTrackerHits", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> _VTXHitCollectionHandle{"VTXTrackerHits", Gaudi::DataHandle::Reader, this};
+
+ DataHandle<edm4hep::TrackCollection> _ClupatraTrackCollectionHandle{"ClupatraTracks", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::TrackCollection> _ClupatraTrackSegmentCollectionHandle{"ClupatraSegmentTracks", Gaudi::DataHandle::Writer, this};
+
+ float _bfield ;
+
+ bool _MSOn ;
+ bool _ElossOn ;
+ bool _SmoothOn ;
+
+
+ int _nRun ;
+ int _nEvt ;
+
+ MarlinTrk::IMarlinTrkSystem* _trksystem ;
+
+ const gear::TPCParameters* _gearTPC ;
+
+ TTree *tree;
+ double omega;
+ double pt;
+ int totalCandidates;
+
+// NNClusterer* _clusterer ;
+
+} ;
+
+#endif
+
+
+
diff --git a/Reconstruction/Tracking/src/Clupatra/NNClusterer.h b/Reconstruction/Tracking/src/Clupatra/NNClusterer.h
new file mode 100644
index 0000000000000000000000000000000000000000..2055d99c6fafbd1c55fc0618527de507300f6e36
--- /dev/null
+++ b/Reconstruction/Tracking/src/Clupatra/NNClusterer.h
@@ -0,0 +1,386 @@
+/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
+#ifndef NNClusterer_h
+#define NNClusterer_h 1
+
+#include <list>
+#include <vector>
+
+// #include "LCRTRelations.h"
+
+/** Nearest neighbour type clusering for arbitrary types.
+ *
+ * @author F.Gaede (DESY)
+ * @version $Id: NNClusterer.h 4052 2012-09-12 09:56:04Z gaede $
+ */
+namespace nnclu {
+
+ /** fast check if integer is in a given range [Min,Max] */
+ template <unsigned Min,unsigned Max >
+ inline bool inRange( int i){ return ( (unsigned int) ( i - Min ) <= (unsigned int) ( Max - Min ) ); }
+
+ /** fast check if integer is not in a given range [Min,Max] */
+ template <unsigned Min,unsigned Max >
+ inline bool notInRange( int i){ return ( (unsigned int) ( i - Min ) > (unsigned int) ( Max - Min ) ); }
+
+
+ // forward declaration:
+ template <class U>
+ class Cluster ;
+
+
+ /** Wrapper class for elements that are clustered, holding a pointer to the actual
+ * object "->first" and a pointer to the cluster this obejct belongs to "->second".
+ *
+ * @see Cluster
+ * @author F.Gaede (DESY)
+ * @version $Id: NNClusterer.h 4052 2012-09-12 09:56:04Z gaede $
+ */
+ template <class T>
+ class Element : public std::pair< T*, Cluster<T>* >{
+
+ typedef T value_type ;
+ typedef Cluster<T> cluster_type ;
+
+ public:
+
+ /** Default c'tor takes a pointer to the original element type object. The optional index can be used to
+ * code nearest neighbour bins, e.g. in z-coordinate to speed up the clustering process.
+ */
+ Element(T* element, int index0 = 0 ) : Index0( index0 ) {
+ Pair::first = element ;
+ Pair::second = 0 ;
+ }
+
+ /** C'tor that also takes a pointer to the cluster this element belongs to - in case seed elements/clusters are used.
+ */
+ Element(T* element , Cluster<T>* cl , int index0 = 0) : Index0( index0 ) {
+ Pair::first = element ;
+ Pair::second = cl ;
+ }
+
+ /** Index that can be used to code nearest neighbour bins, e.g. in z-coordinate
+ * to speed up the clustering process.
+ */
+ int Index0 ;
+
+ protected:
+ typedef std::pair< T*, Cluster<T>* > Pair ;
+
+ /** Don't allow default c'tor w/o element */
+ Element() ;
+ } ;
+
+
+ /** Helper class that creates an Elements for an objects of type T.
+ */
+ template <class T>
+ struct MakeElement{
+ Element<T> operator()( T* t) { return new Element<T>( t ) ; }
+ } ;
+
+
+ /** Extension of std::vector that allows to take ownership of objects pointed to and
+ * delete these when going out of scope.
+ */
+ template <class T>
+ class PtrVector : public std::vector<T*> {
+ typedef std::vector<T*> vec ;
+ bool _isOwner ;
+ public:
+ PtrVector() : _isOwner( false ) {}
+ ~PtrVector() {
+ if( _isOwner )
+ for( typename vec::iterator i = vec::begin(),end = vec::end(); i != end ; delete *i++ ) ; //++i ) delete *i ;
+ }
+ void setOwner( bool val=true ) { _isOwner = val ; }
+ };
+
+
+ /** Extension of std::list that allows to take ownership of objects pointed to and
+ * delete these when going out of scope.
+ */
+ template <class T>
+ class PtrList : public std::list<T*> {
+ typedef std::list<T*> vec ;
+ bool _isOwner ;
+ public:
+ PtrList() : _isOwner( false ) {}
+ ~PtrList() {
+ if( _isOwner )
+ for( typename vec::iterator i = vec::begin(),end = vec::end(); i != end ; delete *i++ ) ; //++i ) delete *i ;
+ }
+ void setOwner( bool val=true ) { _isOwner = val ; }
+ };
+
+
+
+ /** Templated class for generic clusters of Elements that are clustered with
+ * an NN-like clustering algorithm. Effectively this is just a list of elements.
+ *
+ * @see Element
+ * @author F.Gaede (DESY)
+ * @version $Id: NNClusterer.h 4052 2012-09-12 09:56:04Z gaede $
+ */
+ template <class T >
+ class Cluster : public std::list< Element<T> * >{ //, public lcrtrel::LCRTRelations {
+
+ public :
+ typedef Element<T> element_type ;
+ typedef std::list< Element<T> * > base ;
+
+ int ID ; //DEBUG
+
+ Cluster() : ID(0) {}
+
+ /** C'tor that takes the first element */
+ Cluster( Element<T>* element) {
+ static int SID=0 ; //DEBUG
+ ID = SID++ ; //DEBUG
+ addElement( element ) ;
+ }
+
+ /** Add a element to this cluster - updates the element's pointer to cluster */
+ void addElement( Element<T>* element ) {
+
+ element->second = this ;
+ base::push_back( element ) ;
+ }
+
+ // /** Remove all elements from the cluster and reset the cluster association, i.e. elements can be
+ // * used for another clustering procedure.
+ // */
+ // template <class Out>
+ // void takeElements(Out result){
+ // typename Cluster<T>::iterator it = this->begin() ;
+ // while( it != this->end() ){
+ // (*it)->second = 0 ;
+ // result++ = *it ;
+ // it = this->erase(it) ;
+ // }
+ // }
+
+
+ /** Free all elements, ie. reset their cluster association - the elements are still in the list !.
+ */
+ void freeElements(){
+
+ for( typename Cluster<T>::iterator it = this->begin(), end = this->end() ; it != end ; it++ ){
+ (*it)->second = 0 ;
+ }
+
+ // typename Cluster<T>::iterator it = this->begin() ;
+ // while( it != this->end() ){
+ // (*it)->second = 0 ;
+ // it = this->erase(it) ;
+ // }
+ }
+
+ /** Merges all elements from the other cluster cl into this cluster */
+ void mergeClusters( Cluster<T>* cl ) {
+
+ for( typename Cluster<T>::iterator it = cl->begin(), end = cl->end() ; it != end ; it++ ){
+ (*it)->second = this ;
+ }
+ this->merge( *cl ) ;
+ }
+
+ /** D'tor frees all remaining elements that still belong to this cluster */
+ ~Cluster() {
+
+ // typename Cluster<T>::iterator it = this->begin() ;
+ // while( it != this->end() )
+
+ for( typename Cluster<T>::iterator it = this->begin() , end = this->end() ; it != end ; ++it ){
+
+ typename Cluster<T>::value_type h = *it ;
+
+ if( h != 0 && h->second == this )
+ h->second = 0 ;
+
+ }
+ }
+ } ;
+
+
+ template <class T>
+ /** Main class for a nearest neighbour type clustering.
+ *
+ * @author F.Gaede (DESY)
+ * @version $Id: NNClusterer.h 4052 2012-09-12 09:56:04Z gaede $
+ */
+ class NNClusterer{
+
+ public:
+ typedef T value_type ;
+ typedef Cluster<T> cluster_type ;
+ typedef Element<T> element_type ;
+ typedef PtrVector< element_type > element_vector ;
+ typedef PtrVector< cluster_type > cluster_vector ;
+ typedef PtrList< element_type > element_list ;
+ typedef PtrList< cluster_type > cluster_list ;
+
+ /** Simple nearest neighbour (NN) clustering algorithm. Users have to provide an input iterator of
+ * Element objects and an output iterator for the clusters found. The predicate has to have
+ * a method with the following signature: bool operator()( const Element<T>*, const Element<T>*).
+ * All pairs of elements for which this method returns 'true' will be merged into one output cluster
+ * - all other pairs of elements will be in different clusters.
+ */
+
+ template <class In, class Out, class Pred >
+ void cluster( In first, In last, Out result, Pred& pred , const unsigned minSize=1) {
+
+ cluster_vector tmp ;
+ tmp.reserve( 1024 ) ;
+
+ while( first != last ) {
+
+ for( In other = first+1 ; other != last ; other ++ ) {
+
+ if( pred( (*first) , (*other) ) ) {
+
+ if( (*first)->second == 0 && (*other)->second == 0 ) { // no cluster exists
+
+ cluster_type* cl = new cluster_type( (*first) ) ;
+
+ cl->addElement( (*other) ) ;
+
+ tmp.push_back( cl ) ;
+
+ }
+ else if( (*first)->second != 0 && (*other)->second != 0 ) { // two clusters
+
+ if( (*first)->second != (*other)->second ) // don't call merge on identical clusters
+ (*first)->second->mergeClusters( (*other)->second ) ;
+
+ } else { // one cluster exists
+
+ if( (*first)->second != 0 ) {
+
+ (*first)->second->addElement( (*other) ) ;
+
+ } else {
+
+ (*other)->second->addElement( (*first) ) ;
+ }
+ }
+
+ }
+ }
+ ++first ;
+ }
+
+ // remove empty clusters
+ for( typename cluster_vector::iterator i = tmp.begin(); i != tmp.end() ; i++ ){
+
+ if( (*i)->size() > minSize-1 ) {
+
+ result++ = *i ;
+ }
+ else {
+
+ delete *i ;
+ }
+ }
+ }
+
+
+ /** Same as above - but requires the elements to be sorted in index0 (only compare neighbouring bins in index0). */
+
+ template <class In, class Out, class Pred >
+ void cluster_sorted( In first, In last, Out result, Pred& pred , const unsigned minSize=1) {
+
+
+ cluster_vector tmp ;
+ tmp.reserve( 1024 ) ;
+
+ while( first != last ) {
+
+ for( In other = first+1 ; other != last ; other ++ ) {
+
+ // if the elements are sorted we can skip the rest of the inner loop
+ if( notInRange<-1,1>( (*first)->Index0 - (*other)->Index0 ) )
+ break ;
+
+ if( pred( (*first) , (*other) ) ) {
+
+ if( (*first)->second == 0 && (*other)->second == 0 ) { // no cluster exists
+
+ cluster_type* cl = new cluster_type( (*first) ) ;
+
+ cl->addElement( (*other) ) ;
+
+ tmp.push_back( cl ) ;
+
+ }
+ else if( (*first)->second != 0 && (*other)->second != 0 ) { // two clusters
+
+ if( (*first)->second != (*other)->second ) // don't call merge on identical clusters
+ (*first)->second->mergeClusters( (*other)->second ) ;
+
+ } else { // one cluster exists
+
+ if( (*first)->second != 0 ) {
+
+ (*first)->second->addElement( (*other) ) ;
+
+ } else {
+
+ (*other)->second->addElement( (*first) ) ;
+ }
+ }
+
+ }
+ }
+ ++first ;
+ }
+
+ // remove empty clusters
+ for( typename cluster_vector::iterator i = tmp.begin(); i != tmp.end() ; i++ ){
+
+ if( (*i)->size() > minSize-1 ) {
+
+ result++ = *i ;
+ }
+ else {
+
+ delete *i ;
+ }
+ }
+ }
+
+ };
+ //-----------------------------------------------------------------------------------------------------------------------
+
+
+
+ /**Splits a list into two based on a predicate. The new list will
+ * hold all elements for which Pred is true which are in turn removed
+ * from the original list.
+ */
+ template <class List, class Out, class Pred >
+ void split_list( List& list, Out result, Pred pred) {
+
+ typename List::iterator it = list.begin() ;
+
+ while( it != list.end() ){
+
+ if( pred( *it ) ){
+
+ result++ = *it ;
+ it = list.erase( it ) ;
+
+ } else
+ ++it ;
+ }
+ }
+
+
+
+
+ //-----------------------------------------------------------------------------------------------------------------------
+
+} // namespace nnclu
+
+#endif
+
+
diff --git a/Reconstruction/Tracking/src/Clupatra/RuntimeMap.h b/Reconstruction/Tracking/src/Clupatra/RuntimeMap.h
new file mode 100644
index 0000000000000000000000000000000000000000..2b0b14a2ed1b0edcbcd75e427fa942f6d1f19853
--- /dev/null
+++ b/Reconstruction/Tracking/src/Clupatra/RuntimeMap.h
@@ -0,0 +1,16 @@
+#ifndef CLUPATRA_RUNTIMEMAP_H
+#define CLUPATRA_RUNTIMEMAP_H
+#include <map>
+
+template<class U, class V>
+class RuntimeMap {
+ std::map<U, V> data;
+ public:
+ V& operator()(const U& u) {
+ return data[u];
+ }
+ void init() {
+ data.clear();
+ }
+};
+#endif
diff --git a/Reconstruction/Tracking/src/Clupatra/TrackingHelper.h b/Reconstruction/Tracking/src/Clupatra/TrackingHelper.h
new file mode 100644
index 0000000000000000000000000000000000000000..2a16aa8250daccc87c7cd320ea18d90a4e8ed7b1
--- /dev/null
+++ b/Reconstruction/Tracking/src/Clupatra/TrackingHelper.h
@@ -0,0 +1,62 @@
+#ifndef CLUPATRA_TRACKINGHELPER
+#define CLUPATRA_TRACKINGHELPER
+
+#include "edm4hep/TrackState.h"
+#include "edm4hep/Track.h"
+#include "UTIL/BitField64.h"
+#include "UTIL/CellIDDecoder.h"
+#include "UTIL/ILDConf.h"
+#include "lcio.h"
+
+inline bool hasTrackStateAt(edm4hep::ConstTrack track, int location) {
+ for (auto it = track.trackStates_begin(); it != track.trackStates_end(); it++) {
+ if (it->location == location) {
+ return true;
+ }
+ }
+ return false;
+}
+
+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 edm4hep::TrackState();
+}
+
+inline float getTanLambda(const edm4hep::ConstTrack &track) {
+ return track.getTrackStates(0).tanLambda;
+}
+inline float getOmega(const edm4hep::ConstTrack &track) {
+ return track.getTrackStates(0).omega;
+}
+inline float getD0(const edm4hep::ConstTrack &track) {
+ return track.getTrackStates(0).D0;
+}
+inline float getZ0(const edm4hep::ConstTrack &track) {
+ return track.getTrackStates(0).Z0;
+}
+inline float getPhi(const edm4hep::ConstTrack &track) {
+ return track.getTrackStates(0).phi;
+}
+
+
+inline int getLayer(const edm4hep::TrackerHit hit) {
+ UTIL::BitField64* _encoder = new UTIL::BitField64(lcio::ILDCellID0::encoder_string);
+ _encoder->setValue(hit.getCellID());
+ int layer = (*_encoder)[lcio::ILDCellID0::layer];
+ delete _encoder;
+ return layer;
+}
+
+inline int getLayer(const edm4hep::ConstTrackerHit hit) {
+ UTIL::BitField64* _encoder = new UTIL::BitField64(lcio::ILDCellID0::encoder_string);
+ _encoder->setValue(hit.getCellID());
+ int layer = (*_encoder)[lcio::ILDCellID0::layer];
+ delete _encoder;
+ return layer;
+}
+
+#endif
diff --git a/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp b/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..479b981c5803b1112232161cd991eeeaab2194bb
--- /dev/null
+++ b/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
@@ -0,0 +1,1513 @@
+#include "clupatra_new.h"
+#include <set>
+#include <vector>
+
+#include <UTIL/BitField64.h>
+#include <UTIL/ILDConf.h>
+#include <UTIL/BitSet32.h>
+
+///---- GEAR ----
+#include "gear/GEAR.h"
+#include "gear/TPCParameters.h"
+#include "gear/TPCModule.h"
+
+#include "gear/BField.h"
+
+#include "IMPL/TrackerHitImpl.h"
+#include "IMPL/TrackStateImpl.h"
+
+#include "FWCore/DataHandle.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GearSvc/IGearSvc.h"
+using namespace MarlinTrk ;
+
+namespace lcio{
+ // bits 0-15 are reserved !?
+ const int ILDTrackTypeBit::SEGMENT = 16 ;
+ const int ILDTrackTypeBit::COMPOSITE = 17 ;
+}
+
+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;
+
+namespace clupatra_new{
+
+ bool TrackZSort::operator()( edm4hep::ConstTrack l, edm4hep::ConstTrack r){
+ return ( std::abs( TrackInfo_of_edm4hepTrack(l)->zAvg ) < std::abs( TrackInfo_of_edm4hepTrack(r)->zAvg ) ) ;
+ }
+
+ void ComputeTrackerInfo::operator()( edm4hep::ConstTrack o )
+ {
+ edm4hep::ConstTrack lTrk = o ;
+
+ // compute z-extend of this track segment
+ // const edm4hep::TrackerHitVec& hv = lTrk->getTrackerHits() ;
+
+ float zMin = 1e99 ;
+ float zMax = -1e99 ;
+ float zAvg = 0. ;
+
+ if( lTrk.trackerHits_size() > 1 ) {
+ zMin = lTrk.getTrackerHits( 0 ).getPosition()[2] ;
+ zMax = lTrk.getTrackerHits( lTrk.trackerHits_size() -1 ).getPosition()[2] ;
+ zAvg = ( zMax + zMin ) / 2. ;
+ }
+
+ if( zMin > zMax ){ // swap
+ float d = zMax ;
+ zMax = zMin ;
+ zMin = d ;
+ }
+
+ if( ! TrackInfo_of_edm4hepTrack(lTrk) )
+ TrackInfo_of_edm4hepTrack(lTrk) = new TrackInfoStruct ;
+
+ TrackInfo_of_edm4hepTrack(lTrk)->zMin = zMin ;
+ TrackInfo_of_edm4hepTrack(lTrk)->zMax = zMax ;
+ TrackInfo_of_edm4hepTrack(lTrk)->zAvg = zAvg ;
+
+ }
+ bool TrackCircleDistance:: operator()( nnclu::Element<ClupaPlcioConstTrack>* h0, nnclu::Element<ClupaPlcioConstTrack>* h1){
+
+ edm4hep::ConstTrack trk0 = h0->first->edm4hepTrack ;
+ edm4hep::ConstTrack trk1 = h1->first->edm4hepTrack ;
+
+ const TrackInfoStruct* ti0 = TrackInfo_of_edm4hepTrack(trk0);
+ const TrackInfoStruct* ti1 = TrackInfo_of_edm4hepTrack(trk1);
+
+
+ /* FIXME debug Mingrui
+ streamlog_out( DEBUG2 ) << "TrackCircleDistance::operator() : " << trk0->id() << " <-> " << trk1->id()
+ << " ( ti0->zAvg > ti1->zAvg ) = " << ( ti0->zAvg > ti1->zAvg )
+ << std::endl ;
+ */
+
+
+ // don't allow overlaps in z !!!!
+ float epsilon = 0. ;
+ if( ti0->zAvg > ti1->zAvg ){
+
+ if( ti1->zMax > ( ti0->zMin + epsilon ) )
+ return false ;
+
+ } else {
+
+ if( ti0->zMax > ( ti1->zMin + epsilon ) )
+ return false ;
+
+ }
+
+ double tl0 = std::abs( getTanLambda(trk0) ) ;
+ double tl1 = std::abs( getTanLambda(trk1) ) ;
+
+ if( getTanLambda(trk0) * getTanLambda(trk1) < 0. )
+ return false ; // require the same sign
+
+
+ double dtl = 2. * std::abs( tl0 - tl1 ) / ( tl0 + tl1 ) ;
+
+ /* FIXME debug Mingrui
+ streamlog_out( DEBUG2 ) << "TrackCircleDistance::operator() : " << trk0->id() << " <-> " << trk1->id()
+ << " dtl : " << dtl << " std::abs( tl0 + tl1 ) = " << std::abs( tl0 + tl1 )
+ << " ( dtl > 2. * _dCut && std::abs( tl0 + tl1 ) > 1.e-2 ) = " << ( dtl > 2. * _dCut && std::abs( tl0 + tl1 ) > 1.e-2 )
+ << std::endl ;
+ */
+
+ if( dtl > 2. * _dCut && std::abs( tl0 + tl1 ) > 1.e-2 )
+ return false ;
+ // for very steep tracks (tanL < 0.001 ) tanL might differ largely for curlers due to multiple scattering
+
+ double r0 = 1. / getOmega(trk0) ;
+ double r1 = 1. / getOmega(trk1) ;
+
+ double r0abs = std::abs( r0 ) ;
+ double r1abs = std::abs( r1 ) ;
+
+
+ double d0 = getD0(trk0) ;
+ double d1 = getD0(trk1) ;
+
+ // double z0 = trk0->getZ0() ;
+ // double z1 = trk1->getZ0() ;
+
+ double rIP = 20. ;
+
+ // streamlog_out( DEBUG2 ) << "TrackCircleDistance::operator() : " << trk0->id() << " <-> " << trk1->id()
+ // << " ( std::abs( d0 ) < rIP && std::abs( z0 ) < rIP && std::abs( d1 ) < rIP && std::abs( z1 ) < rIP ) "
+ // << ( std::abs( d0 ) < rIP && std::abs( z0 ) < rIP && std::abs( d1 ) < rIP && std::abs( z1 ) < rIP )
+ // << std::endl ;
+ // // // don't merge tracks that come from an area of 20 mm around the IP
+ // if( std::abs( d0 ) < rIP && std::abs( z0 ) < rIP &&
+ // std::abs( d1 ) < rIP && std::abs( z1 ) < rIP )
+ // return false ;
+
+
+ edm4hep::TrackState ts0 = getTrackStateAt(trk0, lcio::TrackState::AtFirstHit ) ;
+ edm4hep::TrackState ts1 = getTrackStateAt(trk1, lcio::TrackState::AtFirstHit ) ;
+
+ /*
+ double z0 = ( ts0 ? ts0->referencePoint()[2] : 0. ) ;
+ double z1 = ( ts1 ? ts1->referencePoint()[2] : 0. ) ;
+ */
+ // FIXME Mingrui since it is no pointer:
+ double z0 = ts0.referencePoint[2];
+ double z1 = ts1.referencePoint[2];
+
+
+ /* FIXME debug Mingrui
+ streamlog_out( DEBUG2 ) << "TrackCircleDistance::operator() : " << trk0->id() << " <-> " << trk1->id()
+ << " ( std::abs( z0 ) < rIP && std::abs( z1 ) < rIP ) "
+ << ( std::abs( z0 ) < rIP && std::abs( z1 ) < rIP )
+ << std::endl ;
+ */
+
+ if( std::abs( z0 ) < rIP && std::abs( z1 ) < rIP )
+ return false ;
+
+
+
+ double p0 = getPhi(trk0) ;
+ double p1 = getPhi(trk1) ;
+
+ double x0 = ( r0 - d0 ) * sin( p0 ) ;
+ double x1 = ( r1 - d1 ) * sin( p1 ) ;
+
+ double y0 = ( d0 - r0 ) * cos( p0 ) ;
+ double y1 = ( d1 - r1 ) * cos( p1 ) ;
+
+ double dr = 2. * std::abs( r0abs - r1abs ) / (r0abs + r1abs ) ;
+
+ double distMS = sqrt ( ( x0 - x1 ) * ( x0 - x1 ) + ( y0 - y1 ) * ( y0 - y1 ) ) ;
+
+ /* FIXME debug Mingrui
+ streamlog_out( DEBUG2 ) << "TrackCircleDistance:: operator() : " << trk0->id() << " <-> " << trk1->id()
+ << "( dr < _dCut * std::abs( r0 ) && distMS < _dCut * std::abs( r0 ) ) "
+ << ( dr < _dCut * std::abs( r0 ) && distMS < _dCut * std::abs( r0 ) )
+ << " dr : " << dr
+ << " _dCut * std::abs( r0 ) " << _dCut * std::abs( r0 )
+ << " distMS :" << distMS
+ << std::endl ;
+ */
+
+ // return ( dr < DRMAX && distMS < _dCut * std::abs( r0 ) ) ;
+ return ( dr < _dCut * std::abs( r0 ) && distMS < _dCut * std::abs( r0 ) ) ;
+
+ }
+
+ /** helper class to compute the chisquared of two points in rho and z coordinate */
+ struct Chi2_RPhi_Z_Hit{
+ // double operator()( const TrackerHit* h, const gear::Vector3D& v1) {
+ double operator()( const ClupaHit* h, const edm4hep::Vector3d& v1_t) {
+
+ gear::Vector3D v1(v1_t[0], v1_t[1], v1_t[2]);
+ // gear::Vector3D v0( h->getPosition()[0] , h->getPosition()[1] , h->getPosition()[2] ) ;
+ const gear::Vector3D& v0 = h->pos ;
+
+ double sigsr = h->edm4hepHit.getCovMatrix()[0] + h->edm4hepHit.getCovMatrix()[2] ;
+ double sigsz = h->edm4hepHit.getCovMatrix()[5] ;
+ // double sigsr = 0.01 ;
+ // double sigsz = 0.1 ;
+
+ double dPhi = std::abs( v0.phi() - v1.phi() ) ;
+ if( dPhi > M_PI )
+ dPhi = 2.* M_PI - dPhi ;
+
+ double dRPhi = dPhi * v0.rho() ;
+
+ double dZ = v0.z() - v1.z() ;
+
+ return dRPhi * dRPhi / sigsr + dZ * dZ / sigsz ;
+ }
+ };
+
+ //-------------------------------------------------------------------------------
+
+ struct CDot{
+ CDot( int i , int j , double d) : I(i) , J(j) , Dot( d ) {}
+ int I ;
+ int J ;
+ double Dot ;
+ };
+
+ bool sort_CDot( const CDot& c0 ,const CDot& c1 ){
+ return c0.Dot > c1.Dot ;
+ }
+
+ //-------------------------------------------------------------------------------
+
+
+ int addHitsAndFilter( CluTrack* clu, HitListVector& hLV , double dChi2Max, double chi2Cut, unsigned maxStep, ZIndex& zIndex, bool backward,
+ MarlinTrk::IMarlinTrkSystem* trkSys ) {
+
+
+ int nHitsAdded = 0 ;
+
+ const double bfield = gearMgr->getBField().at( gear::Vector3D(0.,0.0,0.) ).z() ;
+
+ // Support for more than one module
+ static const gear::TPCParameters* gearTPC = &(gearMgr->getTPCParameters());
+ // The ternary operator is used to make the trick with the static variable which
+ // is supposed to be calculated only once, also for performance reason
+ static const int maxTPCLayerID = (gearMgr->getDetectorName() == "LPTPC" ) ?
+ gearTPC->getModule(0).getNRows() + gearTPC->getModule(2).getNRows() + gearTPC->getModule(5).getNRows() - 1 : // LCTPC
+ gearTPC->getModule(0).getNRows() - 1 ; // ILD
+
+ clu->sort( LayerSortIn() ) ;
+
+ int layer = ( backward ? clu->front()->first->layer : clu->back()->first->layer ) ;
+
+
+ // std::cout << " ====================== addHitsAndFilter(): - layer " << layer << " backward: " << backward << std::endl;
+ // std::cerr << " ====================== addHitsAndFilter(): - layer " << layer << " backward: " << backward << std::endl;
+
+
+ if( layer <= 0 || layer >= maxTPCLayerID )
+ return nHitsAdded ;
+
+
+ Chi2_RPhi_Z_Hit ch2rzh ;
+
+ IMarlinTrack* trk = MarTrkof(clu);
+// {
+// edm4hep::TrackState ts;
+// double chi2; int ndf;
+// trk->getTrackState( ts, chi2, ndf ) ;
+// std::cout << "Omega of this part is: " << ts.omega << " " << chi2 << " " << ndf << std::endl;
+// }
+
+ if( trk == 0 ){
+ // streamlog_out( DEBUG3 ) << " addHitsAndFilter called with null pointer to MarlinTrk - won't do anything " << std::endl ;
+ return nHitsAdded;
+ }
+
+ unsigned step = 0 ;
+
+ UTIL::BitField64 encoder( UTIL::ILDCellID0::encoder_string ) ;
+ encoder[UTIL::ILDCellID0::subdet] = UTIL::ILDDetID::TPC ;
+
+ edm4hep::ConstTrackerHit firstHit; // = 0 ;
+ firstHit.unlink();
+
+ IMarlinTrack* bwTrk = 0 ;
+
+ if( trkSys && backward ) { //==================== only active if called with _trkSystem pointer ============================
+ // std::cout << "It is backward" << std::endl;
+
+ // need to go back in cluster until 4th hit from the start
+ CluTrack::iterator it = clu->begin() , end = clu->end() ;
+ int i=0 ;
+ for( ; it++ != end && i< 3 ; ++i ) ;
+
+ if( it == end ) --it ;
+
+ // std::cout << " ---- addHitsAndFilter : will smooth back to " << i <<"th hit - size of clu " << clu->size() << std::endl ;
+
+ if( ! (*it)->first ){
+
+ // std::cout << " ---- addHitsAndFilter : null pointer at i-th hit : i = " << i << std::endl ;
+ return nHitsAdded ;
+ }
+
+
+ firstHit = (*it)->first->edm4hepHit ;
+
+ int smoothed = trk->smooth( firstHit ) ;
+
+ double chi2 ;
+ int ndf ;
+ edm4hep::TrackState ts ;
+ trk->getTrackState( firstHit , ts, chi2, ndf ) ;
+
+ // std::cout << " -- addHitsAndFilter(): smoothed track segment : " << MarlinTrk::errorCode( smoothed )
+ // << " using track state : " << ts
+ // << " ---- chi2: " << chi2
+ // << " ndf " << ndf
+ // << std::endl ;
+
+
+
+ // std::auto_ptr<MarlinTrk::IMarlinTrack> mTrk( _trksystem->createTrack() ) ;
+ bwTrk = trkSys->createTrack() ;
+
+ //need to add a dummy hit to the track
+ bwTrk->addHit( firstHit ) ; // use the hit we smoothed back to
+
+ // note: backward here is forward in the MarlinTrk sense, ie. in direction of energy loss
+ bwTrk->initialise( ts , bfield , MarlinTrk::IMarlinTrack::forward ) ;
+
+ }
+
+
+ IMarlinTrack* theTrk = ( bwTrk ? bwTrk : trk ) ;
+
+
+ while( step < maxStep + 1 ) {
+
+ layer += ( backward ? +1 : -1 ) ;
+
+ if( layer < 0 || layer > maxTPCLayerID )
+ break ;
+
+
+ encoder[ UTIL::ILDCellID0::layer ] = layer ;
+
+ // gear::Vector3D xv ;
+ edm4hep::Vector3d xv;
+
+ bool hitAdded = false ;
+
+ int layerID = encoder.lowWord() ;
+ int elementID = 0 ;
+
+ // int mode = ( backward ? IMarlinTrack::modeBackward : IMarlinTrack::modeForward ) ;
+ // int mode = ( backward ? IMarlinTrack::modeForward : IMarlinTrack::modeClosest ) ;
+ int mode = IMarlinTrack::modeClosest ;
+
+
+ int intersects = -1 ;
+
+ if( firstHit.isAvailable() ) { // FIXME Mingrui:: Important!! need to check whether firstHit is OK or not.
+// std::cout << "available" << std::endl;
+// std::cout << layerID << std::endl;
+// std::cout << firstHit << std::endl;
+// std::cout << xv << std::endl;
+// std::cout << elementID << std::endl;
+// std::cout << mode << std::endl;
+ intersects = theTrk->intersectionWithLayer( layerID, firstHit, xv, elementID , mode ) ;
+
+ } else {
+
+ intersects = theTrk->intersectionWithLayer( layerID, xv, elementID , mode ) ;
+ }
+
+
+
+ // std::cout << " -- addHitsAndFilter(): looked for intersection - "
+ // << " Step : " << step
+ // << " at layer: " << layer
+ // << " intersects: " << MarlinTrk::errorCode( intersects )
+ // << " next xing point : " << xv ;
+
+
+ if( intersects == IMarlinTrack::success ) { // found a crossing point
+
+ int zIndCP = zIndex.index( xv[2] ) ;
+
+ HitList& hLL = hLV.at( layer ) ;
+
+ double ch2Min = 1.e99 ;
+ Hit* bestHit = 0 ;
+
+ // streamlog_out( DEBUG3 ) << " -- number of hits on layer " << layer << " : " << hLL.size() << std::endl ;
+
+ for( HitList::const_iterator ih = hLL.begin(), end = hLL.end() ; ih != end ; ++ih ){
+
+ // if the z indices differ by more than one we can continue
+ if( nnclu::notInRange<-1,1>( (*ih)->first->zIndex - zIndCP ) )
+ continue ;
+
+ double ch2 = ch2rzh( (*ih)->first , xv ) ;
+
+ if( ch2 < ch2Min ){
+
+ ch2Min = ch2 ;
+ bestHit = (*ih) ;
+ }
+ }//-------------------------------------------------------------------
+
+ if( bestHit != 0 ){
+
+ // streamlog_out( DEBUG3 ) << " ************ bestHit " << bestHit
+ // << " pos : " << (bestHit ? bestHit->first->pos : gear::Vector3D() )
+ // << " chi2: " << ch2Min
+ // << " chi2Cut: " << chi2Cut << std::endl ;
+
+ const gear::Vector3D& hPos = bestHit->first->pos ;
+
+ if( ch2Min < chi2Cut ) {
+
+ double deltaChi = 0. ;
+
+ edm4hep::ConstTrackerHit ht = bestHit->first->edm4hepHit;
+ int addHit = theTrk->addAndFit(ht , deltaChi, dChi2Max ) ;
+
+
+
+ // streamlog_out( DEBUG3 ) << " ***** assigning left over hit : " << errorCode( addHit ) //<< hPos << " <-> " << xv
+ /*
+ << " dist: " << ( hPos - xv ).r()
+ << " chi2: " << ch2Min
+ << " hit errors : rphi=" << sqrt( bestHit->first->edm4hepHit->getCovMatrix()[0]
+ + bestHit->first->edm4hepHit->getCovMatrix()[2] )
+ << " z= " << sqrt( bestHit->first->edm4hepHit->getCovMatrix()[5] )
+ << " ----- deltaChi = " << deltaChi
+ << std::endl ;
+ */
+
+
+ if( addHit == IMarlinTrack::success ) {
+
+ hitAdded = true ;
+
+ hLL.remove( bestHit ) ;
+ clu->addElement( bestHit ) ;
+
+ // firstHit = 0 ; // after we added a hit, the next intersection search should use this last hit...
+ firstHit.unlink();
+
+ ++nHitsAdded ;
+
+ // streamlog_out( DEBUG ) << " ---- track state filtered with new hit ! ------- " << std::endl ;
+ }
+ } // chi2Cut
+ } // bestHit
+
+
+ } else { // no intersection found
+
+ // we stop searching here
+ step = maxStep + 1 ;
+ }
+
+ // reset the step or increment it
+ step = ( hitAdded ? 1 : step + 1 ) ;
+
+ } // while step < maxStep
+
+
+ if( bwTrk ) delete bwTrk ;
+
+ return nHitsAdded ;
+
+ }
+
+ //------------------------------------------------------------------------------------------------------------
+
+ bool addHitAndFilter( int detectorID, int layer, CluTrack* clu, HitListVector& hLV , double dChi2Max, double chi2Cut) {
+
+ Chi2_RPhi_Z_Hit ch2rzh ;
+
+ IMarlinTrack* trk = MarTrkof(clu) ;
+
+ UTIL::BitField64 encoder( UTIL::ILDCellID0::encoder_string ) ;
+
+ encoder[ UTIL::ILDCellID0::subdet ] = detectorID ;
+ encoder[ UTIL::ILDCellID0::layer ] = layer ;
+
+ int layerID = encoder.lowWord() ;
+
+ // gear::Vector3D xv ;
+ edm4hep::Vector3d xv ;
+
+ int elementID = -1 ;
+
+ int intersects = trk->intersectionWithLayer( layerID, xv, elementID , IMarlinTrack::modeClosest ) ;
+
+ // IMarlinTrack::modeBackward , IMarlinTrack::modeForward
+
+ /*
+ streamlog_out( DEBUG2 ) << " ============ addHitAndFilter(): looked for intersection - "
+ << " detector : " << detectorID
+ << " at layer: " << layer
+ << " intersects: " << MarlinTrk::errorCode( intersects )
+ << " next xing point : " << xv ;
+ */
+
+ bool hitAdded = false ;
+
+ if( intersects == IMarlinTrack::success ) { // found a crossing point
+
+ //FIXME: this is copied from ClupatraProcessor
+ static const int ZBins = 160 ;
+ ZIndex zIndex( -2750. , 2750. ,ZBins ) ;
+ int zIndCP = zIndex.index( xv[2] ) ;
+
+ HitList& hLL = hLV.at( layer ) ;
+
+ double ch2Min = 1.e99 ;
+ Hit* bestHit = 0 ;
+
+ for( HitList::const_iterator ih = hLL.begin(), end = hLL.end() ; ih != end ; ++ih ){
+
+ // if the z indices differ by more than one we can continue
+ if( nnclu::notInRange<-1,1>( (*ih)->first->zIndex - zIndCP ) )
+ continue ;
+
+ double ch2 = ch2rzh( (*ih)->first , xv ) ;
+
+ if( ch2 < ch2Min ){
+
+ ch2Min = ch2 ;
+ bestHit = (*ih) ;
+ }
+ }//-------------------------------------------------------------------
+
+
+ // streamlog_out( DEBUG2 ) << " ************ bestHit " << bestHit
+ // << " pos : " << (bestHit ? bestHit->first->pos : gear::Vector3D() )
+ // << " chi2: " << ch2Min
+ // << " chi2Cut: " << chi2Cut << std::endl ;
+
+ if( bestHit != 0 ){
+
+ const gear::Vector3D& hPos = bestHit->first->pos ;
+
+ if( ch2Min < chi2Cut ) {
+
+ double deltaChi = 0. ;
+ edm4hep::ConstTrackerHit bh = bestHit->first->edm4hepHit;
+ int addHit = trk->addAndFit( bh, deltaChi, dChi2Max ) ;
+
+
+
+ // streamlog_out( DEBUG2 ) << " ***** assigning left over hit : " << errorCode( addHit ) //<< hPos << " <-> " << xv
+ /*
+ << " dist: " << ( hPos - xv ).r()
+ << " chi2: " << ch2Min
+ << " hit errors : rphi=" << sqrt( bestHit->first->edm4hepHit->getCovMatrix()[0]
+ + bestHit->first->edm4hepHit->getCovMatrix()[2] )
+ << " z= " << sqrt( bestHit->first->edm4hepHit->getCovMatrix()[5] )
+ << std::endl ;
+ */
+
+
+
+
+
+ if( addHit == IMarlinTrack::success ) {
+
+ hitAdded = true ;
+
+ hLL.remove( bestHit ) ;
+ clu->addElement( bestHit ) ;
+
+
+ // streamlog_out( DEBUG2 ) << " ---- track state filtered with new hit ! ------- " << std::endl ;
+ }
+ } // chi2Cut
+ } // bestHit
+
+ } // intersection found
+
+ return hitAdded ;
+ }
+
+
+ //------------------------------------------------------------------------------------------------------------
+
+ void getHitMultiplicities( CluTrack* clu, std::vector<int>& mult ){
+
+ // int static maxTPCLayerID = marlin::Global::GEAR->getTPCParameters().getPadLayout().getNRows() - 1 ;
+ // HitListVector hLV( maxTPCLayerID ) ;
+ // addToHitListVector( clu->begin(), clu->end(), hLV ) ;
+
+ std::map<int, unsigned > hm ;
+
+ for( CluTrack::iterator it=clu->begin(), end =clu->end() ; it != end ; ++ it ){
+
+ ++hm[ (*it)->first->layer ] ;
+ }
+
+ unsigned maxN = mult.size() - 1 ;
+
+ // for( unsigned i=0 ; i < hLV.size() ; ++i ){
+ // unsigned m = hLV[i].size() ;
+
+ for( std::map<int, unsigned>::iterator it= hm.begin() , end = hm.end() ; it != end ; ++ it ){
+
+ unsigned m = ( it->second < maxN ? it->second : maxN ) ;
+
+ if( m == 0 )
+ continue ;
+
+ ++mult[m] ;
+
+ ++mult[0] ;
+ }
+ }
+
+ //-------------------------------------------------------------------------------------------------------------------------
+ void split_multiplicity( Clusterer::cluster_list& cluList, int layerWithMultiplicity , int N) {
+
+ for( Clusterer::cluster_list::iterator it= cluList.begin(), end= cluList.end() ; it != end ; ++it ){
+
+ CluTrack* clu = *it ;
+
+ std::vector<int> mult( N + 2 ) ;
+ // get hit multiplicities up to N ( N+1 means N+1 or higher )
+ getHitMultiplicities( clu , mult ) ;
+
+
+ // streamlog_out( DEBUG2 ) << " **** split_multiplicity - hit multiplicities: \n" ;
+
+ /*
+ for( unsigned i=0,n=mult.size() ; i<n ; ++i) {
+ streamlog_out( DEBUG2 ) << " m["<<i<<"] = " << mult[i] << "\n" ;
+ }
+ */
+
+
+ if( mult[1] >= layerWithMultiplicity )
+ continue ;
+
+ bool split_cluster = false ;
+
+ for( int m=2 ; m<N ; ++m){
+
+ if( m == 2 && mult[m] >= layerWithMultiplicity ){
+
+ // streamlog_out( DEBUG3 ) << " **** split_multiplicity - create_two_clusters \n" ;
+
+ create_two_clusters( *clu , cluList ) ;
+
+ split_cluster = true ;
+ }
+ else if( m == 3 && mult[m] >= layerWithMultiplicity ){
+
+ // streamlog_out( DEBUG3 ) << " **** split_multiplicity - create_three_clusters \n" ;
+
+ create_three_clusters( *clu , cluList ) ;
+
+ split_cluster = true ;
+ }
+ else if( mult[m] >= layerWithMultiplicity ){
+
+ // streamlog_out( DEBUG3 ) << " **** split_multiplicity - create_n_clusters \n" ;
+
+ create_n_clusters( *clu ,cluList , m ) ;
+
+ split_cluster = true ;
+ }
+
+ }
+ if( split_cluster ){
+
+ clu->clear() ;
+
+ it = cluList.erase( it ) ;
+ --it ; // erase returns iterator to next element
+ }
+
+
+ }
+ }
+
+ //-------------------------------------------------------------------------------------------------------------------------
+
+ void create_n_clusters( Clusterer::cluster_type& hV, Clusterer::cluster_list& cluVec , unsigned n ) {
+
+ if( n < 4 ){
+
+ // streamlog_out( ERROR ) << " create_n_clusters called for n = " << n << std::endl ;
+ return ;
+ }
+
+ // copy of the algorithm for creating three clusters (minimize angle between hits as seen from IP)
+
+ hV.freeElements() ;
+
+ // Support for more than one module
+ //auto _gear = service<IGearSvc>("GearSvc");
+
+ static const gear::TPCParameters* gearTPC = &(gearMgr->getTPCParameters());
+ // The ternary operator is used to make the trick with the static variable which
+ // is supposed to be calculated only once, also for performance reason
+ static const int tpcNRow =
+ (gearMgr->getDetectorName() == "LPTPC" ) ?
+ gearTPC->getModule(0).getNRows() + gearTPC->getModule(2).getNRows() + gearTPC->getModule(5).getNRows() : // LCTPC
+ gearTPC->getModule(0).getNRows() ; // ILD
+
+ HitListVector hitsInLayer( tpcNRow ) ;
+ addToHitListVector( hV.begin(), hV.end(), hitsInLayer ) ;
+
+ std::vector< CluTrack*> clu(n) ;
+
+ std::vector< gear::Vector3D > lastp(n) ;
+
+ for(unsigned i=0; i<n; ++i){
+
+ clu[i] = new CluTrack ;
+
+ cluVec.push_back( clu[i] ) ;
+ }
+
+ for( int l=tpcNRow-1 ; l >= 0 ; --l){
+
+ HitList& hL = hitsInLayer[ l ] ;
+
+ // streamlog_out( DEBUG ) << " create_n_clusters --- layer " << l << " size: " << hL.size() << std::endl ;
+
+ if( hL.size() != n ){ // ignore layers with different hit numbers
+
+ continue ;
+ }
+
+ HitList::iterator iH = hL.begin() ;
+
+ std::vector<Hit*> h(n) ;
+ std::vector< gear::Vector3D> p(n) ;
+
+
+ // streamlog_out( DEBUG2 ) << " create_n_clusters --- layer " << l << std::endl ;
+
+ for(unsigned i=0; i<n; ++i){
+
+ h[i] = *iH++ ;
+
+ p[i] = h[i]->first->pos ;
+
+ // streamlog_out( DEBUG2 ) << " h"<< i << ": " << p[i] << std::endl ;
+ }
+
+ // streamlog_out( DEBUG2 ) << " ---------------------------- " << std::endl ;
+
+
+ if( clu[0]->empty() ){ // first hit tuple
+
+ // streamlog_out( DEBUG ) << " create_n_clusters --- initialize clusters " << std::endl ;
+
+ for(unsigned i=0; i<n; ++i){
+
+ clu[i]->addElement( h[i] ) ;
+
+ lastp[i] = ( 1. / p[i].r() ) * p[i] ;
+ }
+
+ continue ; //---------- that's all for the first layer w/ hits
+ }
+
+
+ // unit vector in direction of current hits
+ std::vector< gear::Vector3D > pu(n) ;
+
+ for(unsigned i=0; i<n; ++i){
+
+ pu[i] = ( 1. / p[i].r() ) * p[i] ;
+ }
+
+ std::list< CDot > cDot ;
+
+ // create list of dot products between last hit in cluster and current hit
+ // cos( angle between hits as seen from IP )
+ for( unsigned i = 0 ; i < n ; ++i ){
+ for( unsigned j = 0 ; j < n ; ++j ){
+
+ cDot.push_back( CDot( i , j , lastp[i].dot( pu[ j ] ) )) ;
+ }
+ }
+
+ // sort dot products in descending order
+ cDot.sort( sort_CDot ) ;
+
+ // assign hits to clusters with largest dot product ( smallest angle )
+
+ std::set<int> cluIdx ;
+ std::set<int> hitIdx ;
+
+ for( std::list<CDot>::iterator it = cDot.begin() ; it != cDot.end() ; ++ it ) {
+ /*
+ streamlog_out( DEBUG2 ) << " I : " << it->I
+ << " J : " << it->J
+ << " d : " << it->Dot
+ << std::endl ;
+ */
+ unsigned i = it->I ;
+ unsigned j = it->J ;
+
+ // ignore clusters or hits that hvae already been assigned
+ if( cluIdx.find( i ) == cluIdx.end() && hitIdx.find( j ) == hitIdx.end() ){
+
+ cluIdx.insert( i ) ;
+ hitIdx.insert( j ) ;
+
+ clu[ i ]->addElement( h[ j ] ) ;
+
+ // cluDir[i] = - ( p[j] - lastp[i] ) ;
+ // cluDir[i] = ( 1. / cluDir[i].r() ) * cluDir[i] ;
+
+ lastp[i] = p[j] ;
+
+ /*
+ streamlog_out( DEBUG2 ) << " **** adding to cluster : " << it->I
+ << " hit : " << it->J
+ << " d : " << it->Dot
+ << std::endl ;
+ */
+ }
+
+ }
+ }
+
+
+ // streamlog_out( DEBUG ) << " create_three_clusters --- clu[0] " << clu[0]->size()
+ // << " clu[1] " << clu[1]->size()
+ // << " clu[2] " << clu[2]->size()
+ // << std::endl ;
+
+
+
+
+ return ;
+ }
+
+
+ //======================================================================================================================
+
+ void create_three_clusters( Clusterer::cluster_type& hV, Clusterer::cluster_list& cluVec ) {
+
+ hV.freeElements() ;
+
+ // auto _gear = service<IGearSvc>("GearSvc");
+ // Support for more than one module
+ static const gear::TPCParameters* gearTPC = &(gearMgr->getTPCParameters());
+ // The ternary operator is used to make the trick with the static variable which
+ // is supposed to be calculated only once, also for performance reason
+ static const int tpcNRow =
+ (gearMgr->getDetectorName() == "LPTPC" ) ?
+ gearTPC->getModule(0).getNRows() + gearTPC->getModule(2).getNRows() + gearTPC->getModule(5).getNRows() : // LCTPC
+ gearTPC->getModule(0).getNRows() ; // ILD
+
+ HitListVector hitsInLayer( tpcNRow ) ;
+ addToHitListVector( hV.begin(), hV.end(), hitsInLayer ) ;
+
+ CluTrack* clu[3] ;
+
+ gear::Vector3D lastp[3] ;
+ // gear::Vector3D cluDir[3] ;
+
+ clu[0] = new CluTrack ;
+ clu[1] = new CluTrack ;
+ clu[2] = new CluTrack ;
+
+ cluVec.push_back( clu[0] ) ;
+ cluVec.push_back( clu[1] ) ;
+ cluVec.push_back( clu[2] ) ;
+
+
+ for( int l=tpcNRow-1 ; l >= 0 ; --l){
+
+ HitList& hL = hitsInLayer[ l ] ;
+
+ // streamlog_out( DEBUG ) << " create_three_clusters --- layer " << l << " size: " << hL.size() << std::endl ;
+
+ if( hL.size() != 3 ){
+
+ continue ;
+ }
+
+ HitList::iterator iH = hL.begin() ;
+
+ Hit* h[3] ;
+ h[0] = *iH++ ;
+ h[1] = *iH++ ;
+ h[2] = *iH ;
+
+ gear::Vector3D p[3] ;
+ p[0] = h[0]->first->pos ;
+ p[1] = h[1]->first->pos ;
+ p[2] = h[2]->first->pos ;
+
+
+ /*
+ streamlog_out( DEBUG ) << " create_three_clusters --- layer " << l
+ << " h0 : " << p[0]
+ << " h1 : " << p[1]
+ << " h2 : " << p[2] ;
+ */
+ // << std::endl ;
+
+ if( clu[0]->empty() ){ // first hit triplet
+
+ // streamlog_out( DEBUG ) << " create_three_clusters --- initialize clusters " << std::endl ;
+
+ clu[0]->addElement( h[0] ) ;
+ clu[1]->addElement( h[1] ) ;
+ clu[2]->addElement( h[2] ) ;
+
+ // lastp[0] = p[0] ;
+ // lastp[1] = p[1] ;
+ // lastp[2] = p[2] ;
+
+ lastp[0] = ( 1. / p[0].r() ) * p[0] ;
+ lastp[1] = ( 1. / p[1].r() ) * p[1] ;
+ lastp[2] = ( 1. / p[2].r() ) * p[2] ;
+
+ continue ;
+ }
+
+
+ // unit vector in direction of current hits
+ gear::Vector3D pu[3] ;
+ pu[0] = ( 1. / p[0].r() ) * p[0] ;
+ pu[1] = ( 1. / p[1].r() ) * p[1] ;
+ pu[2] = ( 1. / p[2].r() ) * p[2] ;
+
+
+ std::list< CDot > cDot ;
+
+ // create list of dot products between last hit in cluster and current hit
+ // cos( angle between hits as seen from IP )
+ for( int i = 0 ; i < 3 ; ++i ){
+ for( int j = 0 ; j < 3 ; ++j ){
+
+ // // unit vector in direction of last hit and
+ // gear::Vector3D pu = - ( p[0] - lastp[0] ) ;
+ // pu[1] = - ( p[1] - lastp[1] ) ;
+ // pu[2] = - ( p[2] - lastp[2] ) ;
+ // pu[0] = ( 1. / pu[0].r() ) * pu[0] ;
+ // pu[1] = ( 1. / pu[1].r() ) * pu[1] ;
+ // pu[2] = ( 1. / pu[2].r() ) * pu[2] ;
+ // cDot.push_back( CDot( i , j , cluDir[i].dot( pu[ j ] ) )) ;
+
+ cDot.push_back( CDot( i , j , lastp[i].dot( pu[ j ] ) )) ;
+
+ }
+ }
+
+ // sort dot products in descending order
+ cDot.sort( sort_CDot ) ;
+
+ // assign hits to clusters with largest dot product ( smallest angle )
+
+ std::set<int> cluIdx ;
+ std::set<int> hitIdx ;
+
+ for( std::list<CDot>::iterator it = cDot.begin() ; it != cDot.end() ; ++ it ) {
+
+ /*
+ streamlog_out( DEBUG ) << " I : " << it->I
+ << " J : " << it->J
+ << " d : " << it->Dot
+ << std::endl ;
+ */
+
+ int i = it->I ;
+ int j = it->J ;
+
+ // ignore clusters or hits that hvae already been assigned
+ if( cluIdx.find( i ) == cluIdx.end() && hitIdx.find( j ) == hitIdx.end() ){
+
+ cluIdx.insert( i ) ;
+ hitIdx.insert( j ) ;
+
+ clu[ i ]->addElement( h[ j ] ) ;
+
+ // cluDir[i] = - ( p[j] - lastp[i] ) ;
+ // cluDir[i] = ( 1. / cluDir[i].r() ) * cluDir[i] ;
+
+ lastp[i] = p[j] ;
+
+ /*
+ streamlog_out( DEBUG ) << " **** adding to cluster : " << it->I
+ << " hit : " << it->J
+ << " d : " << it->Dot
+ << std::endl ;
+ */
+ }
+
+ }
+
+ // lastp[0] = pu[0] ;
+ // lastp[1] = pu[1] ;
+ // lastp[2] = pu[2] ;
+ // lastp[0] = p[0] ;
+ // lastp[1] = p[1] ;
+ // lastp[2] = p[2] ;
+
+
+ }
+
+
+ // streamlog_out( DEBUG ) << " create_three_clusters --- clu[0] " << clu[0]->size()
+ //<< " clu[1] " << clu[1]->size()
+ //<< " clu[2] " << clu[2]->size()
+ //<< std::endl ;
+
+
+
+
+ return ;
+ }
+ //-----------------------------------------------------------------
+
+ void create_two_clusters( Clusterer::cluster_type& clu, Clusterer::cluster_list& cluVec ) {
+
+
+ clu.freeElements() ;
+
+ // streamlog_out( DEBUG ) << " create_two_clusters --- called ! - size : " << clu.size() << std::endl ;
+ // auto _gear = service<IGearSvc>("GearSvc");
+
+
+ // Support for more than one module
+ static const gear::TPCParameters* gearTPC = &(gearMgr->getTPCParameters());
+ // The ternary operator is used to make the trick with the static variable which
+ // is supposed to be calculated only once, also for performance reason
+ static const int tpcNRow =
+ (gearMgr->getDetectorName() == "LPTPC" ) ?
+ gearTPC->getModule(0).getNRows() + gearTPC->getModule(2).getNRows() + gearTPC->getModule(5).getNRows() : // LCTPC
+ gearTPC->getModule(0).getNRows() ; // ILD
+
+ HitListVector hitsInLayer( tpcNRow ) ;
+
+ addToHitListVector( clu.begin(), clu.end(), hitsInLayer ) ;
+
+
+ CluTrack* clu0 = new CluTrack ;
+ CluTrack* clu1 = new CluTrack ;
+
+ cluVec.push_back( clu0 ) ;
+ cluVec.push_back( clu1 ) ;
+
+
+ gear::Vector3D lastDiffVec(0.,0.,0.) ;
+
+ for( int l=tpcNRow-1 ; l >= 0 ; --l){
+
+ HitList& hL = hitsInLayer[ l ] ;
+
+ // streamlog_out( DEBUG ) << " create_two_clusters --- layer " << l << " size: " << hL.size() << std::endl ;
+
+ if( hL.size() != 2 ){
+
+ // if there is not exactly two hits in the layer, we simply copy the unassigned hits to the leftover hit vector
+ // std::copy( hL.begin(), hL.end() , std::back_inserter( hV ) ) ;
+ continue ;
+ }
+
+ HitList::iterator iH = hL.begin() ;
+
+ Hit* h0 = *iH++ ;
+ Hit* h1 = *iH ;
+
+ gear::Vector3D& p0 = h0->first->pos ;
+ gear::Vector3D& p1 = h1->first->pos ;
+
+
+ // streamlog_out( DEBUG ) << " create_two_clusters --- layer " << l
+ //<< " h0 : " << p0
+ // << " h1 : " << p1 ;
+ // << std::endl ;
+
+ if( clu0->empty() ){ // first hit pair
+
+ // streamlog_out( DEBUG ) << " create_two_clusters --- initialize clusters " << std::endl ;
+
+ clu0->addElement( h0 ) ;
+ clu1->addElement( h1 ) ;
+
+ lastDiffVec = p1 - p0 ;
+
+ continue ;
+ }
+
+ gear::Vector3D d = p1 - p0 ;
+
+ float s0 = ( lastDiffVec + d ).r() ;
+ float s1 = ( lastDiffVec - d ).r() ;
+
+ if( s0 > s1 ){ // same orientation, i.e. h0 in this layer belongs to h0 in first layer
+
+ // streamlog_out( DEBUG ) << " create_two_clusters --- same orientation " << std::endl ;
+ clu0->addElement( h0 ) ;
+ clu1->addElement( h1 ) ;
+
+ } else{ // oposite orientation, i.e. h1 in this layer belongs to h0 in first layer
+
+ // streamlog_out( DEBUG2 ) << " create_two_clusters --- oposite orientation " << std::endl ;
+ clu0->addElement( h1 ) ;
+ clu1->addElement( h0 ) ;
+ }
+
+ }
+
+ // clu.freeHits() ;
+
+ // streamlog_out( DEBUG1 ) << " create_two_clusters --- clu0 " << clu0->size()
+ // << " clu1 " << clu1->size() << std::endl ;
+
+ // return std::make_pair( clu0 , clu1 ) ;
+
+ return ;
+ }
+
+ //-------------------------------------------------------------------------------------------------------------------------
+
+
+ // /** Find the nearest hits in the previous and next layers - (at most maxStep layers appart - NOT YET).
+ // */
+ // void findNearestHits( CluTrack& clu, int maxLayerID, int maxStep){
+
+ // HitListVector hitsInLayer( maxLayerID ) ;
+ // addToHitListVector( clu.begin(), clu.end(), hitsInLayer ) ;
+
+ // for(CluTrack::iterator it= clu.begin() ; it != clu.end() ; ++it ){
+
+ // Hit* hit0 = *it ;
+ // gear::Vector3D& pos0 = hit0->first->pos ;
+
+ // int layer = hit0->first->layerID ;
+
+ // int l=0 ;
+ // if( (l=layer+1) < maxLayerID ) {
+
+ // HitList& hL = hitsInLayer[ l ] ;
+
+ // double minDist2 = 1.e99 ;
+ // Hit* bestHit = 0 ;
+
+ // for( HitList::iterator iH = hL.begin() ; iH != hL.end() ; ++iH ) {
+
+ // Hit* hit1 = *iH ;
+ // gear::Vector3D& pos1 = hit1->first->pos ;
+
+ // double dist2 = ( pos0 - pos1 ).r2() ;
+
+ // if( dist2 < minDist2 ){
+ // minDist2 = dist2 ;
+ // bestHit = hit1 ;
+ // }
+ // }
+ // if( bestHit ){
+
+ // hit0->first->nNNHit = bestHit->first ;
+ // hit0->first->nDist = minDist2 ;
+ // }
+ // }
+
+ // if( (l=layer-1) > 0 ) {
+
+ // HitList& hL = hitsInLayer[ l ] ;
+
+ // double minDist2 = 1.e99 ;
+ // Hit* bestHit = 0 ;
+
+ // for( HitList::iterator iH = hL.begin() ; iH != hL.end() ; ++iH ) {
+
+ // Hit* hit1 = *iH ;
+ // gear::Vector3D& pos1 = hit1->first->pos ;
+
+ // double dist2 = ( pos0 - pos1 ).r2() ;
+
+ // if( dist2 < minDist2 ){
+ // minDist2 = dist2 ;
+ // bestHit = hit1 ;
+ // }
+ // }
+ // if( bestHit ){
+
+ // hit0->first->pNNHit = bestHit->first ;
+ // hit0->first->pDist = minDist2 ;
+ // }
+ // }
+
+ // }
+ // }
+
+
+ //-------------------------------------------------------------------------------------------------------------------------
+
+ MarlinTrk::IMarlinTrack* IMarlinTrkFitter::operator() (CluTrack* clu) {
+
+ bool isFirstFit = true ;
+ double maxChi2 = _maxChi2Increment ;
+
+start:
+
+ MarlinTrk::IMarlinTrack* trk = _ts->createTrack();
+
+ //if( clu->empty() ){
+ if( clu->size() < 3 ){
+
+ // streamlog_out( ERROR ) << " IMarlinTrkFitter::operator() : cannot fit cluster track with less than 3 hits ! " << std::endl ;
+
+ return trk ;
+ }
+
+ MarTrkof(clu) = trk ;
+
+ clu->sort( LayerSortOut() ) ;
+
+ // need to reverse the order for incomming track segments (curlers)
+ // assume particle comes from IP
+ Hit* hf = clu->front() ;
+ Hit* hb = clu->back() ;
+
+ bool reverse_order = ( std::abs( hf->first->pos.z() ) > std::abs( hb->first->pos.z()) + 3. ) ;
+
+ unsigned nHit = 0 ;
+
+ if( reverse_order ){
+ // std::cout << "It is true order" << std::endl;
+ for( CluTrack::reverse_iterator it=clu->rbegin() ; it != clu->rend() ; ++it){
+ edm4hep::ConstTrackerHit ph = (*it)->first->edm4hepHit;
+ trk->addHit(ph) ;
+ ++nHit ;
+ // std::cout << " hit added " << (*it)->first->edm4hepHit << std::endl ;
+ }
+
+ trk->initialise( MarlinTrk::IMarlinTrack::forward ) ;
+
+ } else {
+
+ // std::cout << "It is reverse order" << std::endl;
+ for( CluTrack::iterator it=clu->begin() ; it != clu->end() ; ++it){
+ edm4hep::ConstTrackerHit ph = (*it)->first->edm4hepHit;
+ trk->addHit(ph) ;
+ ++nHit ;
+ // std::cout << " hit added "<< (*it)->first->edm4hepHit << std::endl ;
+ }
+
+ trk->initialise( MarlinTrk::IMarlinTrack::backward ) ;
+ }
+
+ int code = trk->fit( maxChi2 ) ;
+ // for (auto hit : hitsInFit) std::cout << hit.first << std::endl;
+
+ if( code != MarlinTrk::IMarlinTrack::success ){
+
+ std::cout << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IMarlinTrkFitter : problem fitting track "
+ << " error code : " << MarlinTrk::errorCode( code )
+ << std::endl ;
+
+ }
+
+
+ std::vector<std::pair<edm4hep::ConstTrackerHit, double> > hitsInFit ;
+ trk->getHitsInFit( hitsInFit ) ;
+ //----- if the fit did not fail but has a small number of hits used,
+ // we try again one more time with a larger max-chi2-increment
+ if( isFirstFit && ( 1.*hitsInFit.size()) / (1.*nHit ) < 0.2 ) { // fixme: parameter
+
+ isFirstFit = false ;
+
+ maxChi2 = 2. * _maxChi2Increment ;
+
+ // streamlog_out( DEBUG4 ) << " >>>>>> IMarlinTrkFitter : small number of hits used in fit " << hitsInFit.size() << "/" << nHit << " = "
+ // << ( 1.*hitsInFit.size()) / (1.*nHit ) << " refit with larger max chi2 increment: " << maxChi2 << std::endl ;
+ delete trk ;
+
+ goto start ; // ;-)
+ }
+ //----------------------------------------------------------------------
+
+
+ return trk;
+ }
+
+
+ //---------------------------------------------------------------------------------------------------------------------------
+
+ edm4hep::ConstTrack PLCIOTrackConverter::operator() (CluTrack* c) {
+
+ static lcio::BitField64 encoder( lcio::ILDCellID0::encoder_string ) ;
+
+ edm4hep::Track trk;
+
+ trk.setType( lcio::ILDDetID::TPC ) ;
+
+ double e = 0.0 ;
+ int nHit = 0 ;
+ for( CluTrack::iterator hi = c->begin(); hi != c->end() ; hi++) {
+
+ // 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 ) ;
+ e += (*hi)->first->edm4hepHit.getEDep() ;
+ nHit++ ;
+ }
+
+ MarlinTrk::IMarlinTrack* mtrk = MarTrkof(c);
+
+ MarTrk_of_edm4hepTrack(trk) = mtrk ;
+
+ trk.setDEdx( ( nHit ? e/nHit : -1. ) ) ;
+
+ /* FIXME Mingrui
+ trk->subdetectorHitNumbers().resize( 2 * lcio::ILDDetID::ETD ) ;
+
+ trk->subdetectorHitNumbers()[ 2*lcio::ILDDetID::TPC - 2 ] = 0 ;
+ trk->subdetectorHitNumbers()[ 2*lcio::ILDDetID::TPC - 1 ] = nHit ;
+ */
+
+ RuntimeMap<edm4hep::ConstTrackerHit, int> DChi2_of_hit;
+
+ if( mtrk != 0 && ! c->empty() ){
+
+
+ std::vector<std::pair<edm4hep::ConstTrackerHit, double> > hitsInFit ;
+ mtrk->getHitsInFit( hitsInFit ) ;
+ // for (auto hit : hitsInFit) std::cout << hit.second << std::endl;
+ // FIXME Mingrui
+ // trk->subdetectorHitNumbers()[ 2*lcio::ILDDetID::TPC - 2 ] = hitsInFit.size() ;
+
+ if( ! hitsInFit.empty() ){
+
+ // store the delta chi2 for the given hit
+ for(unsigned i=0, N=hitsInFit.size() ; i<N ; ++i){
+
+ DChi2_of_hit(hitsInFit[i].first) = hitsInFit[i].second ;
+
+ //IMPL::TrackerHitImpl* thi = dynamic_cast<IMPL::TrackerHitImpl*> ( hitsInFit[i].first ) ;
+ //thi->setQualityBit( UTIL::ILDTrkHitQualityBit::USED_IN_FIT , 1 ) ;
+ }
+
+ edm4hep::TrackState tsIP;
+ edm4hep::TrackState tsFH;
+ edm4hep::TrackState tsLH;
+ edm4hep::TrackState tsCA;
+
+ tsIP.location = lcio::TrackState::AtIP;
+ tsFH.location = lcio::TrackState::AtFirstHit;
+ tsLH.location = lcio::TrackState::AtLastHit;
+ tsCA.location = lcio::TrackState::AtCalorimeter;
+
+ double chi2 ;
+ int ndf ;
+ int code ;
+
+ // Hit* hf = c->front() ;
+ // Hit* hb = c->back() ;
+ // bool reverse_order = ( std::abs( hf->first->pos.z() ) > std::abs( hb->first->pos.z()) + 3. ) ;
+ // lcio::TrackerHit* fHit = ( reverse_order ? hb->first->lcioHit : hf->first->lcioHit ) ;
+ // lcio::TrackerHit* lHit = ( reverse_order ? hf->first->lcioHit : hb->first->lcioHit ) ;
+
+ 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 ========================
+
+ code = mtrk->getTrackState( fHit, tsFH, chi2, ndf ) ;
+
+ if( code != MarlinTrk::IMarlinTrack::success ){
+
+ std::cout << " >>>>>>>>>>> PLCIOTrackConverter : could not get TrackState at first Hit !!?? "
+ << " error code : " << MarlinTrk::errorCode( code )
+ << std::endl ;
+ }
+
+ // ======= get TrackState at last hit ========================
+
+#define use_fit_at_last_hit 0
+
+#if use_fit_at_last_hit
+ 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 );
+// std::cout << "the hit" << std::endl;
+// std::cout << lHit.getCellID0() << std::endl;
+// std::cout << last_constrained_hit << std::endl;
+ 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;
+
+// std::cout << "----------------------------------------------------------" << std::endl;
+// std::cout << "Position" << lHit.getPosition() << std::endl;
+// std::cout << "----------------------------------------------------------" << std::endl;
+
+
+ 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);
+
+// std::cout << "Propagate success" << std::endl;
+
+#endif
+
+// std::cout << "We are here" << std::endl;
+ if( code != MarlinTrk::IMarlinTrack::success ){
+
+ std::cout << " >>>>>>>>>>> PLCIOTrackConverter : could not get TrackState at last Hit !!?? " << std::endl ;
+ }
+
+ // ======= get TrackState at calo face ========================
+ //
+ encoder.reset() ;
+ encoder[ lcio::ILDCellID0::subdet ] = lcio::ILDDetID::ECAL ;
+ encoder[ lcio::ILDCellID0::layer ] = 0 ;
+ encoder[ lcio::ILDCellID0::side ] = lcio::ILDDetID::barrel;
+ int layerID = encoder.lowWord() ;
+ int sensorID = -1 ;
+
+#if use_fit_at_last_hit
+ code = mtrk->propagateToLayer( layerID , lHit, tsCA, chi2, ndf, sensorID, MarlinTrk::IMarlinTrack::modeClosest ) ;
+#else // get the track state at the calorimter from a propagating from the last(first) constrained fit position
+ code = mtrk->propagateToLayer( layerID , last_constrained_hit, tsCA, chi2, ndf, sensorID, MarlinTrk::IMarlinTrack::modeClosest ) ;
+#endif
+
+// std::cout << "We are here 2" << std::endl;
+ if( code == MarlinTrk::IMarlinTrack::no_intersection ){
+
+ encoder[ lcio::ILDCellID0::side ] = ( lHit.getPosition()[2] > 0. ? lcio::ILDDetID::fwd : lcio::ILDDetID::bwd ) ;
+ layerID = encoder.lowWord() ;
+
+#if use_fit_at_last_hit
+ code = mtrk->propagateToLayer( layerID , lHit, tsCA, chi2, ndf, sensorID, MarlinTrk::IMarlinTrack::modeClosest ) ;
+#else // get the track state at the calorimter from a propagating from the last(first) constrained fit position
+ code = mtrk->propagateToLayer( layerID , last_constrained_hit, tsCA, chi2, ndf, sensorID, MarlinTrk::IMarlinTrack::modeClosest ) ;
+#endif
+
+ }
+ if ( code !=MarlinTrk::IMarlinTrack::success ) {
+
+ // streamlog_out( DEBUG6 ) << " >>>>>>>>>>> PLCIOTrackConverter : could not get TrackState at calo face !!?? " << std::endl ;
+ }
+
+ // ======= get TrackState at IP ========================
+
+// std::cout << "We are here 3" << std::endl;
+ // const gear::Vector3D ipv( 0.,0.,0. );
+ const double tmp_ipv[] = {0, 0, 0};
+ const edm4hep::Vector3d ipv(tmp_ipv);
+
+ // fg: propagate is quite slow and might not really be needed for the TPC
+
+ code = ( UsePropagate ? mtrk->propagate( ipv, fHit, tsIP, chi2, ndf ) : mtrk->extrapolate( ipv, tsIP, chi2, ndf ) ) ;
+
+ if( code != MarlinTrk::IMarlinTrack::success ){
+
+ std::cout << " >>>>>>>>>>> PLCIOTrackConverter : could not extrapolate TrackState to IP !!?? " << std::endl ;
+ }
+
+ trk.addToTrackStates( tsIP ) ;
+ trk.addToTrackStates( tsFH ) ;
+ trk.addToTrackStates( tsLH ) ;
+ trk.addToTrackStates( tsCA ) ;
+
+// std::cout << "We are here 3" << std::endl;
+ double RMin = sqrt( tsFH.referencePoint[0] * tsFH.referencePoint[0]
+ + tsFH.referencePoint[1] * tsFH.referencePoint[1] ) ;
+
+ trk.setRadiusOfInnermostHit( RMin ) ;
+
+ trk.setChi2( chi2 ) ;
+ trk.setNdf( ndf ) ;
+
+ } else {
+
+ std::cout << " >>>>>>>>>>> PLCIOTrackConverter::operator() - hitsInFitEmpty ! - nHits " << nHit << std::endl ;
+ }
+
+ } else {
+
+ // this is not an error for debug collections that just consist of track hits (no fit )
+ // streamlog_out( ERROR ) << " >>>>>>>>>>> LCIOTrackConverter::operator() (CluTrack* c) : "
+ // << " no MarlinTrk::IMarlinTrack* found for cluster !!?? " << std::endl ;
+ }
+
+
+ return trk;
+ }
+
+
+ }//namespace
diff --git a/Reconstruction/Tracking/src/Clupatra/clupatra_new.h b/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
new file mode 100644
index 0000000000000000000000000000000000000000..5e5c9529ec48547d4ea960c7c811314664340017
--- /dev/null
+++ b/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
@@ -0,0 +1,598 @@
+#ifndef clupatra_new_h
+#define clupatra_new_h
+
+#include "RuntimeMap.h"
+#include "TrackingHelper.h"
+
+#include <cmath>
+#include <algorithm>
+#include <time.h>
+#include <math.h>
+#include <sstream>
+#include <memory>
+#include "assert.h"
+
+#include "NNClusterer.h"
+
+#include "lcio.h"
+#include "EVENT/TrackerHit.h"
+#include "IMPL/TrackImpl.h"
+#include "UTIL/Operators.h"
+#include "UTIL/CellIDDecoder.h"
+#include "UTIL/ILDConf.h"
+
+#include "gear/GEAR.h"
+
+#include "edm4hep/TrackState.h"
+
+#include "TrackSystemSvc/IMarlinTrack.h"
+#include "TrackSystemSvc/IMarlinTrkSystem.h"
+
+// ----- include for verbosity dependend logging ---------
+// #include "marlin/VerbosityLevels.h"
+
+
+/** Helper structs that should go to LCIo to make extraction of layer (and subdetector etc easier )
+*/
+namespace lcio{
+ struct ILDDecoder : public CellIDDecoder<TrackerHit>{
+ ILDDecoder() : lcio::CellIDDecoder<TrackerHit>( lcio::ILDCellID0::encoder_string ) {}
+ } ;
+ static ILDDecoder ILD_cellID ;
+
+ struct ILDTrackTypeBit{
+ static const int SEGMENT ;
+ static const int COMPOSITE ;
+ } ;
+}
+
+namespace clupatra_new{
+
+ /** Small wrapper extension of the LCIO Hit
+ */
+ struct ClupaHit {
+
+ ClupaHit() :layer(-1),
+ zIndex(-1),
+ phiIndex(-1),
+ edm4hepHit(0),
+ pos(0.,0.,0.) {}
+ int layer ;
+ int zIndex ;
+ int phiIndex ;
+ edm4hep::ConstTrackerHit edm4hepHit ;
+ gear::Vector3D pos ;
+
+ };
+
+ // inline lcio::TrackerHit* lcioHit( const ClupaHit* h) { return h->lcioHit ; }
+
+
+ //------------------ typedefs for elements and clusters ---------
+
+ typedef nnclu::NNClusterer< ClupaHit > Clusterer ;
+
+ typedef Clusterer::element_type Hit ;
+ typedef Clusterer::cluster_type CluTrack ;
+
+ typedef Clusterer::element_vector HitVec ;
+ typedef Clusterer::cluster_vector CluTrackVec ;
+
+ typedef std::list<Hit*> HitList ;
+ typedef std::vector< HitList > HitListVector ;
+
+ struct ClupaPlcioConstTrack {
+ edm4hep::ConstTrack edm4hepTrack ;
+ ClupaPlcioConstTrack(edm4hep::ConstTrack edm4hepTrack) : edm4hepTrack(edm4hepTrack) {}
+ };
+
+ // typedef GenericHitVec<ClupaHit> GHitVec ;
+ // typedef GenericClusterVec<ClupaHit> GClusterVec ;
+
+ //------------------------------------------------------------------------------------------
+
+ // struct GHit : lcrtrel::LCExtension<GHit, Hit > {} ;
+
+ //------------------------------------------------------------------------------------------
+
+ // struct MarTrk : lcrtrel::LCExtension<MarTrk, MarlinTrk::IMarlinTrack> {} ;
+
+ //------------------------------------------------------------------------------------------
+
+ // struct DChi2 : lcrtrel::LCFloatExtension<DChi2> {} ;
+
+ //----------------------------------------------------------------
+
+ /** Simple predicate class for computing an index from N bins of the z-coordinate of LCObjects
+ * that have a float/double* getPostion() method.
+ */
+ class ZIndex{
+ public:
+ /** C'tor takes zmin and zmax - NB index can be negative and larger than N */
+ ZIndex( float zmin , float zmax , int n ) : _zmin( zmin ), _zmax( zmax ) , _N (n) {}
+
+ template <class T>
+ inline int operator() (T* hit) {
+
+ return (int) std::floor( ( hit->getPosition()[2] - _zmin ) / ( _zmax - _zmin ) * _N ) ;
+ }
+
+ inline int index( double z) { return (int) std::floor( ( z - _zmin ) / ( _zmax - _zmin ) * _N ) ; }
+
+ protected:
+ ZIndex() {} ;
+ float _zmin ;
+ float _zmax ;
+ int _N ;
+ } ;
+
+ //------------------------------------------------------------------------------------------
+
+ struct ZSort {
+ inline bool operator()( const Hit* l, const Hit* r) {
+ return ( l->first->pos.z() < r->first->pos.z() );
+ }
+ };
+
+
+ //------------------------------------------------------------------------------------------
+
+ struct PtSort { // sort tracks wtr to pt - largest first
+ inline bool operator()( const edm4hep::ConstTrack l, const edm4hep::ConstTrack r) {
+ return ( std::abs( getOmega(l) ) < std::abs( getOmega(r) ) ); // pt ~ 1./omega
+ }
+ };
+
+ //------------------------------------------------------------------------------------------
+
+
+ /** Add the elements (Hits) from (First,Last) to the HitListVector - vector needs to be initialized, e.g. with
+ * hLV.resize( MAX_LAYER_INDEX )
+ */
+ template <class First, class Last>
+ void addToHitListVector( First first, Last last, HitListVector& hLV ){
+
+ while( first != last ){
+
+ // FIXME debug Mingrui
+ //streamlog_out( DEBUG ) << " add hit << " << *first << " to layer " << (*first)->first->layer << std::endl ;
+
+ hLV[ (*first)->first->layer ].push_back( *first ) ;
+ ++first ;
+ }
+ }
+
+ //------------------------------------------------------------------------------------------
+
+ /** Predicate class for 'distance' of NN clustering. */
+ class HitDistance{
+ public:
+
+ HitDistance(float dCut, float caCut = -1.0 ) : _dCutSquared( dCut*dCut ) , _caCut( caCut ) {}
+
+ /** Merge condition: true if distance is less than dCut */
+ inline bool operator()( Hit* h0, Hit* h1){
+
+ if( std::abs( h0->Index0 - h1->Index0 ) > 1 ) return false ;
+
+ if( h0->first->layer == h1->first->layer )
+ return false ;
+
+ if( _caCut > 0. && std::abs( h0->first->layer - h1->first->layer ) == 1 ){
+
+ gear::Vector3D& p0 = h0->first->pos ;
+ gear::Vector3D& p1 = h1->first->pos ;
+
+ double cosAlpha = p0.dot( p1 ) / p0.r() / p1.r() ;
+
+ // merge hits that seem to come from stiff track from the IP
+ //fixme: make parameter
+ if( cosAlpha > _caCut ) return true ;
+ }
+
+ return ( h0->first->pos - h1->first->pos).r2() < _dCutSquared ;
+ }
+ protected:
+ HitDistance() ;
+ float _dCutSquared, _caCut ;
+ } ;
+
+
+ // /** Predicate class for 'distance' of NN clustering. */
+
+ // struct HitDistance{ float _dCutSquared ;
+ // HitDistance(float dCut) : _dCutSquared( dCut*dCut ) {}
+
+ // inline bool operator()( Hit* h0, Hit* h1){
+
+ // if( h0->first->layer == h1->first->layer )
+ // return false ;
+
+ // return ( h0->first->pos - h1->first->pos).r2()
+ // < _dCutSquared ;
+ // }
+ // } ;
+
+
+
+ //-----------------------------------------------
+
+ struct PLCIOTrackConverter{
+
+ bool UsePropagate ;
+ PLCIOTrackConverter() : UsePropagate(false ) {}
+
+ edm4hep::ConstTrack operator() (CluTrack* c) ;
+
+ } ;
+
+ //------------------------------------------------------------------------------------------
+ /** Predicate class for identifying small clusters. */
+ struct ClusterSize {
+ ClusterSize( unsigned n) : _n(n) {}
+ bool operator()(const CluTrack* cl) const { return cl->size() < _n ; }
+ unsigned _n ;
+ };
+
+ //------------------------------------------------------------------------------------------
+ /** Predicate class for identifying clusters with duplicate pad rows - returns true
+ * if the fraction of duplicate hits is larger than 'fraction'.
+ */
+ struct DuplicatePadRows{
+
+ unsigned _N ;
+ float _f ;
+ DuplicatePadRows(unsigned nPadRows, float fraction) : _N( nPadRows), _f( fraction ) {}
+
+ bool operator()(const CluTrack* cl) const {
+
+ // check for duplicate layer numbers
+ std::vector<int> hLayer( _N ) ;
+ typedef CluTrack::const_iterator IT ;
+
+ unsigned nHit = 0 ;
+ for(IT it=cl->begin() ; it != cl->end() ; ++it ) {
+ ++ hLayer[ (*it)->first->layer] ;
+ ++ nHit ;
+ }
+ unsigned nDuplicate = 0 ;
+ for(unsigned i=0 ; i < _N ; ++i ) {
+ if( hLayer[i] > 1 )
+ nDuplicate += hLayer[i] ;
+ }
+ return double(nDuplicate)/nHit > _f ;
+ }
+ };
+
+ //------------------------------------------------------------------------------------------
+
+ // helper for sorting cluster wrt layerID
+ struct LayerSortOut{
+ bool operator()( const Hit* l, const Hit* r) { return l->first->layer < r->first->layer ; }
+ } ;
+ struct LayerSortIn{
+ bool operator()( const Hit* l, const Hit* r) { return l->first->layer > r->first->layer ; }
+ } ;
+
+
+
+
+ //------------------------------------------------------------------------------------------
+
+ struct IMarlinTrkFitter{
+
+ MarlinTrk::IMarlinTrkSystem* _ts ;
+ double _maxChi2Increment ;
+
+ IMarlinTrkFitter(MarlinTrk::IMarlinTrkSystem* ts, double maxChi2Increment=DBL_MAX ) :
+ _ts( ts ) ,
+ _maxChi2Increment(maxChi2Increment) {}
+
+
+ MarlinTrk::IMarlinTrack* operator() (CluTrack* clu) ;
+ };
+
+ //-------------------------------------------------------------------------------------
+
+ /** Try to add hits from hLV (hit lists per layer) to the cluster. The cluster needs to have a fitted KalTrack associated to it.
+ * Hits are added if the resulting delta Chi2 is less than dChiMax - a maxStep is the maximum number of steps (layers) w/o
+ * successfully merging a hit.
+ */
+ int addHitsAndFilter( CluTrack* clu, HitListVector& hLV , double dChiMax, double chi2Cut, unsigned maxStep, ZIndex& zIndex, bool backward=false,
+ MarlinTrk::IMarlinTrkSystem* trkSys=0) ;
+ //------------------------------------------------------------------------------------------
+
+ /** Try to add a hit from the given HitList in layer of subdetector to the track.
+ * A hit is added if the resulting delta Chi2 is less than dChiMax.
+ */
+ bool addHitAndFilter( int detectorID, int layer, CluTrack* clu, HitListVector& hLV , double dChiMax, double chi2Cut) ;
+
+ //------------------------------------------------------------------------------------------
+ /** Split up clusters that have a hit multiplicity of 2,3,4,...,N in at least layersWithMultiplicity.
+ */
+ void split_multiplicity( Clusterer::cluster_list& cluList, int layersWithMultiplicity , int N=5) ;
+
+ //------------------------------------------------------------------------------------------
+ /** Returns the number of rows where cluster clu has i hits in mult[i] for i=1,2,3,4,.... -
+ * mult[0] counts all rows that have hits
+ */
+ void getHitMultiplicities( CluTrack* clu, std::vector<int>& mult ) ;
+
+ //------------------------------------------------------------------------------------------
+
+ /** Split the cluster into two clusters.
+ */
+ void create_two_clusters( Clusterer::cluster_type& clu, Clusterer::cluster_list& cluVec ) ;
+
+ //------------------------------------------------------------------------------------------
+
+ /** Split the cluster into three clusters.
+ */
+ void create_three_clusters( Clusterer::cluster_type& clu, Clusterer::cluster_list& cluVec ) ;
+
+
+ /** Split the cluster into N clusters.
+ */
+ void create_n_clusters( Clusterer::cluster_type& clu, Clusterer::cluster_list& cluVec , unsigned n ) ;
+
+
+
+ //=======================================================================================
+
+ struct TrackInfoStruct{
+ TrackInfoStruct() : zMin(0.), zAvg(0.), zMax(0.), startsInner(false), isCentral(false), isForward(false), isCurler(false) {}
+ float zMin ;
+ float zAvg ;
+ float zMax ;
+ bool startsInner ;
+ bool isCentral ;
+ bool isForward ;
+ bool isCurler ;
+ } ;
+ // struct TrackInfo : lcrtrel::LCOwnedExtension<TrackInfo, TrackInfoStruct> {} ;
+ typedef TrackInfoStruct TrackInfo;
+ /** Helper class to compute track segment properties.
+ */
+ struct ComputeTrackerInfo{
+ void operator()( edm4hep::ConstTrack o );
+ };
+
+ //=======================================================================================
+
+ /** helper class for merging split track segments */
+
+ class TrackSegmentMerger{
+
+ public:
+ /** C'tor takes merge distance */
+ TrackSegmentMerger(float chi2Max, MarlinTrk::IMarlinTrkSystem* trksystem, float b) : _chi2Max( chi2Max ) , _trksystem( trksystem), _b(b) {}
+
+ float _chi2Max ;
+ MarlinTrk::IMarlinTrkSystem* _trksystem ;
+ float _b ;
+
+ /** Merge condition: ... */
+ inline bool operator()( nnclu::Element<ClupaPlcioConstTrack>* h0, nnclu::Element<ClupaPlcioConstTrack>* h1){
+
+ edm4hep::ConstTrack trk0 = h0->first->edm4hepTrack ;
+ edm4hep::ConstTrack trk1 = h1->first->edm4hepTrack ;
+
+
+ // protect against merging multiple segments (and thus complete tracks)
+ if( h0->second || h1->second )
+ return false ;
+
+ // const TrackInfoStruct* ti0 = trk0->ext<TrackInfo>() ;
+ // const TrackInfoStruct* ti1 = trk1->ext<TrackInfo>() ;
+
+ // const lcio::TrackState* tsF0 = trk0->getTrackState( lcio::TrackState::AtFirstHit ) ;
+ // const lcio::TrackState* tsL0 = trk0->getTrackState( lcio::TrackState::AtLastHit ) ;
+
+ // const lcio::TrackState* tsF1 = trk1->getTrackState( lcio::TrackState::AtFirstHit ) ;
+ // const lcio::TrackState* tsL1 = trk1->getTrackState( lcio::TrackState::AtLastHit ) ;
+
+ // --- get three hits per track : first last midddle
+
+ unsigned nhit0 = trk0.trackerHits_size() ;
+ unsigned nhit1 = trk1.trackerHits_size() ;
+
+ edm4hep::ConstTrackerHit thf0 = trk0.getTrackerHits( 0 ) ;
+ edm4hep::ConstTrackerHit thf1 = trk1.getTrackerHits( 0 ) ;
+
+ edm4hep::ConstTrackerHit thl0 = trk0.getTrackerHits( nhit0 - 1 ) ;
+ edm4hep::ConstTrackerHit thl1 = trk1.getTrackerHits( nhit1 - 1 ) ;
+
+ // lcio::TrackerHit* thm1 = trk1->getTrackerHits()[ nhit1 / 2 ] ;
+ // lcio::TrackerHit* thm0 = trk0->getTrackerHits()[ nhit0 / 2 ] ;
+
+ // int lthf0 = ILD_cellID( thf0 )[ lcio::ILDCellID0::layer ] ;
+ // int lthf1 = ILD_cellID( thf1 )[ lcio::ILDCellID0::layer ] ;
+
+ // int lthl0 = ILD_cellID( thl0 )[ lcio::ILDCellID0::layer ] ;
+ // int lthl1 = ILD_cellID( thl1 )[ lcio::ILDCellID0::layer ] ;
+ int lthf0 = getLayer(thf0);
+ int lthf1 = getLayer(thf1);
+ int lthl0 = getLayer(thl0);
+ int lthl1 = getLayer(thl1);
+
+ // if( lthf0 <= lthl1 && lthf1 <= lthl0 ) return false ;
+
+ // allow the track segements to overlap slightly - FIXME: make a parameter ...
+ const int overlapRows = 4 ;
+ if( lthf0 + overlapRows <= lthl1 && lthf1 + overlapRows <= lthl0 ) return false ;
+
+ // 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 ) ;
+
+ bool outward = ( nhit0 > nhit1 ? lthl0 <= lthf1 + overlapRows : lthl1 <= lthf0 + overlapRows ) ;
+
+ unsigned n = oth.trackerHits_size() ;
+
+ edm4hep::ConstTrackerHit th0 = ( outward ? oth.getTrackerHits(0) : oth.getTrackerHits(n-1) ) ;
+ edm4hep::ConstTrackerHit th1 = (oth.getTrackerHits(n/2) );
+ edm4hep::ConstTrackerHit th2 = ( outward ? oth.getTrackerHits(n-1) : oth.getTrackerHits(0) );
+
+
+ // track state at last hit migyt be rubish....
+ // const lcio::TrackState* ts = ( outward ? trk->getTrackState( lcio::TrackState::AtLastHit ) : trk->getTrackState( lcio::TrackState::AtFirstHit ) ) ;
+ const edm4hep::TrackState ts = getTrackStateAt(trk, lcio::TrackState::AtFirstHit ) ;
+
+
+
+ // FIXME Mingrui debug
+ // streamlog_out( DEBUG3 ) << " ******* TrackSegmentMerger : will extrapolate track " << ( outward ? " outwards\t" : " inwards\t" )
+ //<< lcio::lcshort( trk ) << " vs: [" << std::hex << oth->id() << std::dec << "]" << std::endl ;
+
+ // if( trk->id() == 0x0004534d && oth->id() == 0x000454a6 ){
+ // streamlog_out( DEBUG3 ) << " &&&&&&&&&&&&&& Track 1 : \n" << *trk
+ // << " &&&&&&&&&&&&&& Track 2 : \n" << *oth
+ // << std::endl ;
+ // }
+
+ std::auto_ptr<MarlinTrk::IMarlinTrack> mTrk( _trksystem->createTrack() ) ;
+
+ int nHit = trk.trackerHits_size() ;
+
+ if( nHit == 0 /*|| ts ==0*/ ) // FIXME Mingrui Since it is no pointer
+ return false ;
+
+ // float initial_chi2 = trk->getChi2() ;
+ // float initial_ndf = trk->getNdf() ;
+
+ // FIXME Mingrui debug
+ // streamlog_out( DEBUG3 ) << " -- extrapolate TrackState : " << lcshort( ts ) << std::endl ;
+
+ edm4hep::ConstTrackerHit ht = trk.getTrackerHits(0);
+ //need to add a dummy hit to the track
+ mTrk->addHit( ht ) ; // is this the right hit ??????????
+
+ mTrk->initialise( ts , _b , MarlinTrk::IMarlinTrack::backward ) ;
+
+ double deltaChi ;
+ int addHit = 0 ;
+
+ //----- now try to add the three hits : ----------------
+ addHit = mTrk->addAndFit( th0 , deltaChi, _chi2Max ) ;
+
+ /*
+ * FIXME Mingrui debug
+ streamlog_out( DEBUG3 ) << " **** adding first hit : " << gear::Vector3D( th0->getPosition() )
+ << " added : " << MarlinTrk::errorCode( addHit )
+ << " deltaChi2: " << deltaChi
+ << std::endl ;
+ */
+
+ if( addHit != MarlinTrk::IMarlinTrack::success ) return false ;
+
+ //---------------------
+ addHit = mTrk->addAndFit( th1 , deltaChi, _chi2Max ) ;
+
+ /* FIXME Mingrui debug
+ streamlog_out( DEBUG3 ) << " **** adding second hit : " << gear::Vector3D( th1->getPosition() )
+ << " added : " << MarlinTrk::errorCode( addHit )
+ << " deltaChi2: " << deltaChi
+ << std::endl ;
+ */
+
+ if( addHit != MarlinTrk::IMarlinTrack::success ) return false ;
+
+ //--------------------
+ addHit = mTrk->addAndFit( th2 , deltaChi, _chi2Max ) ;
+
+ /* FIXME Mingrui debug
+ streamlog_out( DEBUG3 ) << " **** adding third hit : " << gear::Vector3D( th2->getPosition() )
+ << " added : " << MarlinTrk::errorCode( addHit )
+ << " deltaChi2: " << deltaChi
+ << std::endl ;
+ */
+
+ if( addHit != MarlinTrk::IMarlinTrack::success ) return false ;
+
+
+
+ ////////////
+ return true ;
+ ///////////
+ }
+
+ };
+ //=======================================================================================
+
+ /** helper class for merging track segments, based on circle (and tan lambda) */
+
+ class TrackCircleDistance{
+
+ public:
+ /** C'tor takes merge distance */
+ TrackCircleDistance(float dCut) : _dCutSquared( dCut*dCut ) , _dCut(dCut){}
+
+ /** Merge condition: ... */
+ bool operator() ( nnclu::Element<ClupaPlcioConstTrack>* h0, nnclu::Element<ClupaPlcioConstTrack>* h1);
+ protected:
+ float _dCutSquared ;
+ float _dCut ;
+ //=======================================================================================
+ };
+
+ struct TrackZSort { // sort tracks wtr to abs(z_average )
+ bool operator()( edm4hep::ConstTrack l, edm4hep::ConstTrack r);
+ };
+
+
+ //=======================================================================================
+
+ /** Helper class that creates an Elements for an LCOjects of type T.
+ */
+ template <class T>
+ struct MakePLCIOElement{
+ nnclu::Element<T>* operator()( T* o) { return new nnclu::Element<T>(o) ; }
+ } ;
+
+
+ //=======================================================================================
+
+ class Timer{
+ public:
+ Timer(){
+ _clocks.reserve( 100 ) ;
+ _names.reserve( 100 ) ;
+
+ _clocks.push_back(0) ;
+ _names.push_back( "start" ) ;
+ }
+ unsigned registerTimer( const std::string& name ){
+ _clocks.push_back(0) ;
+ _names.push_back( name ) ;
+ return _clocks.size() - 1 ;
+ }
+
+ void time(unsigned index){
+ _clocks[ index ] = clock() ;
+ }
+ void start() { time(0) ; }
+
+
+ std::string toString(){
+
+ std::stringstream s ;
+
+ s << " ============= Timer ================================ " << std::endl ;
+ unsigned N=_clocks.size() ;
+ for( unsigned i=1 ; i < N ; ++i){
+ s << " " << _names[i] << " : " << double(_clocks[i] - _clocks[i-1] ) / double(CLOCKS_PER_SEC) << std::endl ;
+ }
+ s << " Total : " << double(_clocks[N-1] - _clocks[0] ) / double(CLOCKS_PER_SEC) << std::endl ;
+ s << " ==================================================== " << std::endl ;
+
+ return s.str() ;
+ }
+ protected:
+ std::vector< clock_t> _clocks ;
+ std::vector< std::string > _names ;
+ };
+
+
+}
+#endif