.. _program_listing_file_larflow_CRTMatch_CRTHitMatch.cxx:

Program Listing for File CRTHitMatch.cxx
========================================

|exhale_lsh| :ref:`Return to documentation for file <file_larflow_CRTMatch_CRTHitMatch.cxx>` (``larflow/CRTMatch/CRTHitMatch.cxx``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include "CRTHitMatch.h"
   
   #include "LArUtil/LArProperties.h"
   #include "larflow/Reco/cluster_functions.h"
   
   #include "TRandom2.h"
   
   namespace larflow {
   namespace crtmatch {
   
     void CRTHitMatch::clear() {
   
       // input containers
       _intime_opflash_v.clear();
       _outtime_opflash_v.clear();
       _crthit_v.clear();
       _crttrack_v.clear();
       _lfcluster_v.clear();
       _pcaxis_v.clear();
       _hit2track_rank_v.clear();
       _all_rank_v.clear();
   
       // output containers
       _matched_hitidx.clear();
       _flash_v.clear();
       _match_hit_v.clear();
       _matched_cluster.clear();
       _matched_cluster_t.clear();    
       _matched_opflash_v.clear();
       _used_tracks_v.clear();
       
     }
   
     void CRTHitMatch::addIntimeOpFlashes( const larlite::event_opflash& opflash_v ) {
       for ( auto const& opf : opflash_v )
         _intime_opflash_v.push_back( opf );
     }
   
     void CRTHitMatch::addCosmicOpFlashes( const larlite::event_opflash& opflash_v ) {
       for ( auto const& opf : opflash_v )
         _outtime_opflash_v.push_back( opf );
     }
   
     void CRTHitMatch::addCRThits( const larlite::event_crthit& crthit_v ) {
       for ( auto const& hit : crthit_v )
         _crthit_v.push_back( hit );
     }
   
     void CRTHitMatch::addCRTtracks( const larlite::event_crttrack& crttrack_v ) {
       for ( auto const& track : _crttrack_v )
         _crttrack_v.push_back( track );
     }
   
     void CRTHitMatch::addLArFlowClusters( const larlite::event_larflowcluster& lfcluster_v,
                                           const larlite::event_pcaxis& pcaxis_v ) {
       
       for ( auto const& cluster : lfcluster_v ) {
         _lfcluster_v.push_back( &cluster );
       }
       
       for ( auto const& pca : pcaxis_v ) {
         _pcaxis_v.push_back( &pca );
       }
       
     }
   
     bool CRTHitMatch::process( larcv::IOManager& iolcv, larlite::storage_manager& llio ) {
   
       clear();
   
       larlite::event_opflash* beamflash
         = (larlite::event_opflash*)llio.get_data( larlite::data::kOpFlash, "simpleFlashBeam" );
       larlite::event_opflash* cosmicflash
         = (larlite::event_opflash*)llio.get_data( larlite::data::kOpFlash, "simpleFlashCosmic" );
       addIntimeOpFlashes( *beamflash );
       addCosmicOpFlashes( *cosmicflash );
   
       // get crt hits
       larlite::event_crthit* crthit_v
         = (larlite::event_crthit*)llio.get_data( larlite::data::kCRTHit, "crthitcorr" );
       addCRThits( *crthit_v );
   
       // get clusters
       larlite::event_larflowcluster* lfclusters_v
         = (larlite::event_larflowcluster*)llio.get_data( larlite::data::kLArFlowCluster,
                                                          _input_cluster_treename );
       larlite::event_pcaxis* pcaxis_v
         = (larlite::event_pcaxis*)llio.get_data( larlite::data::kPCAxis,
                                                  _input_pcaxis_treename );
       addLArFlowClusters( *lfclusters_v, *pcaxis_v );
   
       std::vector< larflow::reco::cluster_t > test_v;
       int icluster = 0;
       for ( auto const& plfcluster : _lfcluster_v ) {
         std::cout << "[convert cluster " << icluster  << "]" << std::endl;
         larflow::reco::cluster_t c = larflow::reco::cluster_from_larflowcluster( *plfcluster );
         test_v.emplace_back( std::move(c) );
         icluster++;
       }
       LARCV_INFO() << "passed conversion test" << std::endl;
       
       return makeMatches();
     }
   
     bool CRTHitMatch::makeMatches() {
   
       // compile matches
       compilematches();
   
       _used_tracks_v.resize( _lfcluster_v.size(), 0 );
   
       // process matches in greedy fashion
       for ( auto& m : _all_rank_v ) {
   
         // get index of crt hit of this match
         int crthitidx = m.hitidx;
   
         auto& hit_matches_v = _hit2track_rank_v[ crthitidx ];
   
         if ( hit_matches_v.size()==0 ) continue; // shouldn't happen
         
         if ( hit_matches_v.size()==1 ) {
   
           // if track already used. skip
           if ( _used_tracks_v[ m.trackidx ]==1 )
             continue;
   
           // if not, make the match
           LARCV_INFO() << " store single crt-hit to cluster match products" << std::endl;
           _matched_hitidx.push_back( crthitidx );
           // create a copy
           larlite::larflowcluster  lfc   = *_lfcluster_v[ m.trackidx ]; // create a copy
           _matched_cluster.emplace_back( std::move(lfc) );
           _used_tracks_v[ m.trackidx ] = 1;
   
         }
         else {
           // more than 1, attempt a merge
   
           // set the usage vector for these matches
           bool has_unused = false;
           std::vector<int> used_track_v( hit_matches_v.size(), 0 );
           for ( size_t i=0; i<hit_matches_v.size(); i++ ) {
             used_track_v[i] = _used_tracks_v[ hit_matches_v[i].trackidx ];
             if ( used_track_v[i]==0 )
               has_unused = true;
           }
   
           if ( has_unused ) {
             bool performed_merge = false;          
             larlite::larflowcluster merge = _merge_matched_cluster( hit_matches_v,
                                                                     used_track_v,
                                                                     performed_merge );
   
             // make match based on merge result
             if (performed_merge) {
               larflow::reco::cluster_t clout = larflow::reco::cluster_from_larflowcluster( merge );            
               _matched_hitidx.push_back( crthitidx );
               _matched_cluster.push_back( merge );
               _used_tracks_v[ m.trackidx ] = 1;          
               LARCV_INFO() << "store merged clusters" << std::endl;
             }
           }
   
         }
         
         
       }//end of loop over all rank matches
   
       LARCV_NORMAL() << "number of crt-hit to track matches: " << _matched_cluster.size() << std::endl;
   
       // apply t0 offset to clusters
       for ( size_t imatch=0; imatch<_matched_hitidx.size(); imatch++ ) {
         int hitidx = _matched_hitidx[ imatch ];
         auto const& crthit = _crthit_v[ hitidx ];
         auto&  lfc = _matched_cluster[ imatch ];
         float xoffset = larutil::LArProperties::GetME()->DriftVelocity()*( crthit.ts2_ns*0.001 );
         // we need to t0 subtract the points
         for ( auto& lfhit : lfc ) {
           lfhit[0] -= xoffset;
         }
         // store cluster_t for modified larflow cluster
         larflow::reco::cluster_t clout = larflow::reco::cluster_from_larflowcluster( lfc );
         _matched_cluster_t.emplace_back( std::move(clout) );      
       }
   
       // now match CRT hits with tracks to opflashes
       // make list of flashes
       for ( auto const& flash : _intime_opflash_v )
         _flash_v.push_back( &flash );
       for ( auto const& flash : _outtime_opflash_v )
         _flash_v.push_back( &flash );
       // make list of crt hits with matches
       for ( auto const& hitidx : _matched_hitidx ) {
         _match_hit_v.push_back( &_crthit_v[ hitidx ] );
       }
       _matchOpflashes( _flash_v, _match_hit_v, _matched_cluster, _matched_opflash_v );
   
       return true;
       
     }
   
     void CRTHitMatch::compilematches() {
   
       _hit2track_rank_v.resize( _crthit_v.size() );
   
       // plane positions
       float crt_plane_pos[4] = { -261.606, -142.484, 393.016, 618.25 }; // bottom, side, side, top
   
       //larlite::event_pcaxis* ev_pcaout = (larlite::event_pcaxis*)outio.get_data( larlite::data::kPCAxis, "crtmatch" );
   
       // loop over larflow tracks, assign it to the best hit in each CRT plane
       // we are filling out _hit2track_rank_v[ hit index ] -> list of track indices
       LARCV_INFO() << "==============================" << std::endl;
       LARCV_INFO() << " Start match function" << std::endl;
       if ( _pcaxis_v.size()!=_lfcluster_v.size() ) {
         throw std::runtime_error( "number of clusters and pc axis do not match" );
       }
       
       for (size_t itrack=0; itrack<_pcaxis_v.size(); itrack++ ) {
   
         // represent track and its direction with its pc axis
         auto const& pca = *_pcaxis_v[itrack];
   
         // enforce minimum length
         float len = getLength(pca);
         if ( len<5.0 ) continue;
   
         // allocate vars for best intersection for each CRT plane
         float min_dist[4]  = { 1.0e9, 1.0e9, 1.0e9, 1.0e9 };
         int best_hitidx[4] = {-1, -1, -1, -1};
         std::vector<float> best_panel_pos[4];
         for (int i=0; i<4; i++) best_panel_pos[i].resize(3,0);
   
         // loop over crt hits
         for (size_t jhit=0; jhit<_crthit_v.size(); jhit++) {
           auto const& crthit = _crthit_v[jhit];
           int crtplane = crthit.plane;
   
           // calculate distance to crt hit and position of intersection on the panel
           std::vector<float> panel_pos;
           float dist = makeOneMatch( pca, crthit, panel_pos );
   
           // update closest hit on the plane
           if ( dist>0 && dist<min_dist[crtplane] ) {
             min_dist[crtplane] = dist;
             best_hitidx[crtplane] = jhit;
             best_panel_pos[crtplane] = panel_pos;
           }
         }//end of hit loop
         LARCV_DEBUG() << " [" << itrack << "] closest dist per plane = "
                       << "[ " << min_dist[0] << ", " << min_dist[1] << ", " << min_dist[2] << ", " << min_dist[3] << "]"
                       << std::endl;
   
         // loop over plane
         for (int p=0; p<4; p++ ) {
   
           // if a good hit found in the plane
           if ( best_hitidx[p]>=0 ) {
             auto const& besthit = _crthit_v[best_hitidx[p]];        
             LARCV_DEBUG() << " crt_hit=(" << besthit.x_pos << ", " << besthit.y_pos << "," << besthit.z_pos << ") "
                           << " panel_pos=(" << best_panel_pos[p][0] << "," << best_panel_pos[p][1] << "," << best_panel_pos[p][2] << ") "
                           << std::endl;
   
             if ( min_dist[p]>=0 && min_dist[p]<50.0 ) {
               float line[3] = {0};
               float dist = 0.;
               for (int i=0; i<3; i++ ) {
                 line[i] = best_panel_pos[p][i]-pca.getAvePosition()[i];
                 dist += line[i]*line[i];
               }
               dist = sqrt(dist);
               for (int i=0; i<3; i++ ) line[i] /= dist;
               
               // store a match
               match_t m;
               m.hitidx   = best_hitidx[p];
               m.trackidx = itrack;
               m.tracklen = len;
               m.dist2hit = min_dist[p];
               
               _hit2track_rank_v[ best_hitidx[p] ].push_back( m );
               _all_rank_v.push_back( m );
               
               // store for visualization purposes
               // larlite::pcaxis::EigenVectors e_v; // 3-axis + 2-endpoints
               // for ( size_t p=0; p<3; p++ ) {
               //   std::vector<double> da_v = { line[0], line[1], line[2] };
               //   e_v.push_back( da_v );
               // }
               // std::vector<double> centroid_v = { pca.getAvePosition()[0], pca.getAvePosition()[1], pca.getAvePosition()[2] };
               // std::vector<double> panel_v = { best_panel_pos[p][0], best_panel_pos[p][1], best_panel_pos[p][2] };
               // e_v.push_back( centroid_v );
               // e_v.push_back( panel_v );
               // double eigenval[3] = { min_dist[p], 0, 0 };
               // larlite::pcaxis llpca( true, 1, eigenval, e_v, centroid_v.data(), 0, itrack );
               // //ev_pcaout->emplace_back( std::move(llpca) );
               
             } // close enough to a hit to register a match
           }//if valid hit found for plane
         }//end of loop over CRT planes
         
       }//end of loop over tracks
   
       // sort the all rank list
       std::sort( _all_rank_v.begin(), _all_rank_v.end() );
   
       // sort each hit list
       for ( auto& hitmatches : _hit2track_rank_v ) {
         std::sort( hitmatches.begin(), hitmatches.end() );
       }
   
       printHitInfo();
       
       return;
     }
   
     void CRTHitMatch::printHitInfo() {
       
       LARCV_NORMAL() << "===============================================" << std::endl;
       LARCV_NORMAL() << "[ CRT HIT INFO ]" << std::endl;
       for ( size_t idx=0; idx<_crthit_v.size(); idx++ ) {
         auto const& hit = _crthit_v.at(idx);
         std::stringstream ss;
         ss << " [hit " << idx
            << " p=" << hit.plane
            << " pos=(" << hit.x_pos << ", " << hit.y_pos << ", " << hit.z_pos << ") "
            << " t=" << hit.ts2_ns*1.0e-3 << " usec ]";
         
         ss << " matches[";
         for ( auto const& hidx : _hit2track_rank_v[idx] )
           ss << " (" << hidx.dist2hit << "," << hidx.tracklen << ")";
         ss << " ] ";
         LARCV_NORMAL() << ss.str() << std::endl;
       }
       LARCV_NORMAL() << "===============================================" << std::endl;    
       
     }
   
     float CRTHitMatch::makeOneMatch( const larlite::pcaxis& lfcluster_axis,
                                      const larlite::crthit& hit,
                                      std::vector<float>& panel_pos ) {
       
       int   crt_plane_dim[4] = {        1,        0,       0,      1 }; // Y-axis, X-axis, X-axis, Y-axis
       float endpts[2][3];
       float center[3];
       float dir[3];
       float crtpos[3] = { hit.x_pos, hit.y_pos, hit.z_pos };
       float len = 0;
       for ( int i=0; i<3; i++ ) {
         endpts[0][i] = lfcluster_axis.getEigenVectors()[3][i];
         endpts[1][i] = lfcluster_axis.getEigenVectors()[4][i];
         center[i]    = lfcluster_axis.getAvePosition()[i];
         dir[i] = endpts[0][i]-center[i];
         len += dir[i]*dir[i];
       }
       if (len>0) {
         len = sqrt(len);
         for (int i=0; i<3; i++) dir[i] /= len;
       }
   
       float max_x = 0.;
       float min_x = 1.0e9;
       for (int i=0; i<2; i++ ) {
         if ( endpts[i][0] > max_x )
           max_x = endpts[i][0];
         if ( endpts[i][0] < min_x )
           min_x = endpts[i][0];
       }
       
       // check if parallel to the plane
       if ( dir[ crt_plane_dim[ hit.plane ] ]==0 ) {
         return -1;
       }
   
       // check if in time
       float x_offset = hit.ts2_ns*0.001*larutil::LArProperties::GetME()->DriftVelocity();
       if ( max_x-x_offset >= 260.0 || max_x-x_offset<-10.0 ) {
         return -1;
       }
       if ( min_x-x_offset >= 260.0 || min_x-x_offset<-10.0 ) {
         return -1;
       }
   
       // in time, so let's calculte the position onto the plane
       // remove the time offset
       center[0] -= x_offset;
   
       // project onto crt plane
       float s = (crtpos[ crt_plane_dim[ hit.plane ] ] - center[ crt_plane_dim[hit.plane] ])/dir[ crt_plane_dim[hit.plane] ];
       panel_pos.resize(3,0);
       for ( int i=0; i<3; i++ ) {
         panel_pos[i] = center[i] + s*dir[i];
       }
   
       // calculate distance to the hit
       float dist = 0.;
       for (int i=0; i<3; i++) {
         dist += ( crtpos[i]-panel_pos[i] )*( crtpos[i]-panel_pos[i] );
       }
       dist = sqrt( dist );
   
       panel_pos[0] += x_offset;
   
       return dist;
     }
   
     float CRTHitMatch::getLength( const larlite::pcaxis& pca ) {
       float dist = 0.;
       for (int i=0; i<3; i++ ) {
         float dx = ( pca.getEigenVectors()[3][i]-pca.getEigenVectors()[4][i] );
         dist += dx*dx;
       }
       return sqrt(dist);
     }
   
     void CRTHitMatch::_matchOpflashes( const std::vector< const larlite::opflash* >& flash_v,
                                        const std::vector< const larlite::crthit* >&  hit_v,
                                        const std::vector< larlite::larflowcluster >& cluster_v,
                                        std::vector< larlite::opflash >& matched_opflash_v ) {
   
       std::vector< int > used_flash_v( flash_v.size(), 0 );
       const float _max_dt_flash_crt = 2.0;
   
       for ( size_t ihit=0; ihit<hit_v.size(); ihit++ ) {
   
         auto const& crt = *hit_v[ihit];
           
         std::vector< const larlite::opflash* > matched_in_time_v;
         std::vector<float> dt_usec_v;
         std::vector< int > matched_index_v;
   
         float closest_time = 10e9;
   
         for ( size_t i=0; i<flash_v.size(); i++ ) {
           if ( used_flash_v[i]>0 ) continue;
             
           float dt_usec = crt.ts2_ns*1.0e-3 - flash_v[i]->Time();
   
           if ( fabs(dt_usec) < closest_time )
             closest_time = fabs(dt_usec);
   
           if ( fabs(dt_usec)<_max_dt_flash_crt ) {
             matched_in_time_v.push_back( flash_v[i] );
             dt_usec_v.push_back( dt_usec );
             matched_index_v.push_back( i );
             LARCV_INFO() << "  ...  crt-track and opflash matched. dt_usec=" << dt_usec << std::endl;
           }
             
         }//end of flash loop
   
         LARCV_NORMAL() << "crt-hit has " << matched_index_v.size() << " flash matches. closest time=" << closest_time  << std::endl;
   
         
         if ( matched_in_time_v.size()==0 ) {
           // make empty opflash at the time of the crttrack
           std::vector< double > PEperOpDet(32,0.0);
           larlite::opflash blank_flash( crt.ts2_ns*1.0e-3,
                                         0.0,
                                         crt.ts2_ns*1.0e-3,
                                         0,
                                         PEperOpDet );
           matched_opflash_v.emplace_back( std::move(blank_flash) );
         }
         else if ( matched_in_time_v.size()==1 ) {
           // store a copy of the single matched opflash
           matched_opflash_v.push_back( *matched_in_time_v[0] );
           used_flash_v[ matched_index_v[0] ] = 1;
         }
         else {
   
           // choose closest using distance from center of flash and track
           // need mean of cluster
           std::vector<float> cluster_mean(3,0);
           for ( auto const& hit : cluster_v[ ihit ] ) {
             for (int i=0; i<3; i++)
               cluster_mean[i] += hit[i];
           }
           for (int i=0; i<3; i++ )
             cluster_mean[i] /= (float)cluster_v[ihit].size();
           
           float smallest_dist = 1.0e9;
           const larlite::opflash* closest_opflash = nullptr;
   
           // we have a loose initial standard. if more than one, we use a tighter standard
           bool have_tight_match = false;
           for ( auto& dt_usec : dt_usec_v ) {
             if ( fabs(dt_usec)<1.5 ) have_tight_match = true;
           }
   
           for ( size_t i=0; i<matched_in_time_v.size(); i++ ) {
   
             // ignore loose match if we know we have at least one good one
             if ( have_tight_match && dt_usec_v[i]>1.5 ) continue;
   
             // get mean of opflash
             std::vector<float> flashcenter = { 0.0,
                                                (float)matched_in_time_v[i]->YCenter(),
                                                (float)matched_in_time_v[i]->ZCenter() };
   
             float dist = 0.;
             for (int v=1; v<3; v++ ) {
               dist += ( flashcenter[v]-cluster_mean[v] )*( flashcenter[v]-cluster_mean[v] );
             }
             dist = sqrt(dist);
   
             if ( dist<smallest_dist ) {
               smallest_dist = dist;
               closest_opflash = matched_in_time_v[i];
             }
             LARCV_INFO() << "  ... distance between opflash and track center: " << dist << " cm" << std::endl;
           }//end of candidate opflash match
           
           // store best match
           if ( closest_opflash ) {
             matched_opflash_v.push_back( *closest_opflash );
           }
           else {
             // shouldnt get here
             throw std::runtime_error( "should not get here" );
           }
         }//else more than 1 match
             
       }//end of CRT HIT loop
   
       // check we have the right amount of flashes
       if ( matched_opflash_v.size()!=cluster_v.size() || matched_opflash_v.size()!=hit_v.size() ) {
         throw std::runtime_error( "different amount of input crt-hit, cluster, and opflashes");
       }
       
     }
   
     larlite::larflowcluster CRTHitMatch::_merge_matched_cluster( const std::vector< CRTHitMatch::match_t >& hit_match_v,
                                                                  std::vector<int>& used_in_merge,
                                                                  bool& merged ) {
   
       merged = false;
       larlite::larflowcluster lfcluster;
   
       // if one or zero matches, nothing to do here.
       if ( hit_match_v.size()<=1 )
         return lfcluster;
   
       // we use cluster tools, so start by converting best track 
       larflow::reco::cluster_t mergecluster;
       int nmerged = 0;
       for ( size_t idx=0; idx<hit_match_v.size(); idx++ ) {
         LARCV_INFO() << "merger: make cluster for larflowcluster index=" << hit_match_v[idx].trackidx << " nhits=" << _lfcluster_v[ hit_match_v[idx].trackidx ]->size() << std::endl;
         larflow::reco::cluster_t c = larflow::reco::cluster_from_larflowcluster( *_lfcluster_v[ hit_match_v[idx].trackidx ] );
         if ( nmerged==0 ) {
           // first cluster
           LARCV_INFO() << "merger: set first cluster" << std::endl;
           std::swap(mergecluster,c);
           used_in_merge[idx] = 1;
           merged = true;
           nmerged++;
           continue;
         }
   
         // determine a good match
         LARCV_INFO() << "merger: test for merger" << std::endl;      
         std::vector< std::vector<float> > closestpts_vv;
         float endptdist = larflow::reco::cluster_closest_endpt_dist( mergecluster, c, closestpts_vv );
         float cospca = cluster_cospca( mergecluster, c );
         cospca = 1.0 - fabs(cospca);
   
         if ( float(cospca)*180.0/3.14159 < 20.0 && endptdist < 30.0 ) {
           LARCV_INFO() << "merge test passed" << std::endl;      
           larflow::reco::cluster_t testmerge = larflow::reco::cluster_merge( mergecluster, c );
   
           if ( testmerge.pca_eigenvalues[1] < 30.0 ) {
             merged = true;
             used_in_merge[idx] = 1;
             nmerged++;
             LARCV_INFO() << "swap for merged cluster" << std::endl;                
             std::swap( mergecluster, testmerge );
           }
           
         }
       }
   
       for ( size_t idx=0; idx<hit_match_v.size(); idx++ ) {
         if ( used_in_merge[idx]==1 ) {
           for ( auto const& hit : *_lfcluster_v[ hit_match_v[idx].trackidx ] )
             lfcluster.push_back( hit );
         }
       }
   
       return lfcluster;
     }
   
     void CRTHitMatch::save_to_file( larlite::storage_manager& outll, bool remove_if_no_flash ) {
   
       // larlite::event_crthit* out_crthit
       //   = (larlite::event_crthit*)outll.get_data( larlite::data::kCRTHit, "matchcrthit" );
       larlite::event_crttrack* out_crttrack
         = (larlite::event_crttrack*)outll.get_data( larlite::data::kCRTTrack, "matchcrthit" );
       larlite::event_opflash* out_opflash
         = (larlite::event_opflash*)outll.get_data( larlite::data::kOpFlash, "matchcrthit" );
       larlite::event_larflowcluster* out_lfcluster
         = (larlite::event_larflowcluster*)outll.get_data( larlite::data::kLArFlowCluster, "matchcrthit" );
   
       for (size_t i=0; i<_matched_cluster.size(); i++ ) {
         auto& cluster = _matched_cluster[i];
         auto& opflash = _matched_opflash_v[i];
         auto const& crthit  = *_match_hit_v[i];
         auto const& clustert = _matched_cluster_t[i];
         
         if ( remove_if_no_flash && opflash.TotalPE()==0.0 ) {
           LARCV_INFO() << "no matching flash for matched CRT hit[" << i << "], not saving" << std::endl;
           continue;
         }
         
         float petot = opflash.TotalPE();
         LARCV_NORMAL() << "saving track with opflash pe=" << petot << " nopdets=" << opflash.nOpDets() << std::endl;
   
         // form a crt-track
         larlite::crttrack outtrack;
         // outtrack.feb_id = crthit.feb_id;
         // outtrack.pesmap = crthit.pesmap;
         // outtrack.peshit = crthit.peshit;
         // outtrack.ts0_s       = crthit.ts0_s;
         // outtrack.ts0_ns      = crthit.ts0_ns;
         // outtrack.ts0_ns_corr = crthit.ts0_ns_corr;
   
         // first crt hit
         outtrack.ts2_ns_h1  = crthit.ts2_ns;
         outtrack.plane1     = crthit.plane;
         outtrack.x1_pos     = crthit.x_pos;
         outtrack.y1_pos     = crthit.y_pos;
         outtrack.z1_pos     = crthit.z_pos;
   
         // second crt hit: use the farthest point on the cluster
         outtrack.plane2     = -1;
         outtrack.ts2_ns_h2 = crthit.ts2_ns;
   
         std::vector<float> crtpos = { outtrack.x1_pos, outtrack.y1_pos, outtrack.z1_pos };
   
         int farend = 0;
         float fardist = 0.;
         for ( int iend=0; iend<2; iend++ ) {
           float dist = 0;
           for (int i=0; i<3; i++ ) {
             dist += ( crtpos[i]- clustert.pca_ends_v[iend][i] )*( crtpos[i]- clustert.pca_ends_v[iend][i] );
           }
           dist = sqrt(dist);
           if ( dist>fardist ) {
             fardist = dist;
             farend = iend;
           }
         }
   
         outtrack.x2_pos = clustert.pca_ends_v[farend][0];
         outtrack.y2_pos = clustert.pca_ends_v[farend][1];
         outtrack.z2_pos = clustert.pca_ends_v[farend][2];      
         
         out_crttrack->push_back(outtrack);
         out_opflash->push_back(opflash);
         out_lfcluster->push_back( cluster );
         
       }
       
     }
   
     bool CRTHitMatch::was_cluster_used( int idx ) {
       if (idx<0 || idx>=_used_tracks_v.size() )
         return false;
   
       if ( _used_tracks_v[idx]==1 )
         return true;
   
       return false;
     }
   
   }
   }