---

mMuiRoadFinder1.cxx

Go to the documentation of this file.

#include <iostream>
#include <cmath>

#include "PHNode.h"
#include "PHCompositeNode.h"
#include "PHDataNode.h"
#include "PHIODataNode.h"
#include "PHNodeIterator.h"

#include "mMuiRoadFinder1.h"
#include "utiLineFit.h"
#include "TMuiRoadOGroup.h"
#include "MuiGeomClasses.hh"

#include <list>
#include <algorithm>
#include <vector>

//MUIOO interfaces
#include <TMutNode.h>

struct muioo_cluster_distance
{
  TMuiClusterMapO::pointer pCluster;
  float distance;
};

class muioo_cluster_distance_sort{  //sort in order of distance
public:
  bool operator() (const muioo_cluster_distance&, const muioo_cluster_distance&) const;
};

bool muioo_cluster_distance_sort::operator()(const muioo_cluster_distance& a, const muioo_cluster_distance& b) const
{
  if(a.distance<b.distance) return true;
  if(a.distance>b.distance) return false;
  if(a.distance==b.distance){
    if(a.distance<b.distance) return true;
  }
  return false;
}


// Default constructor and destructor to pacify CINT
mMuiRoadFinder1::mMuiRoadFinder1() :
   _roadmap(0),
   _road1dmap(0),
   _clustermap(0),
   _timer( PHTimeServer::get()->insert_new("mMuiRoadFinder1") )

{
  name = "mMuiRoadFinder1";
}

mMuiRoadFinder1::~mMuiRoadFinder1(){}

PHBoolean
mMuiRoadFinder1::event(PHCompositeNode *top_node) {

  _timer.get()->restart();

#ifdef MUIOO_DEBUG
  cout<<"Begin mMuiRoadFinder1."<<endl;
#endif
  set_interface_ptrs(top_node);

  find_1droads();

  cut_1droads();

  find_2droads();

  flag_golden();
  
#ifdef MUIOO_DEBUG  
  TMui1DRoadMapO::iterator iroad = _road1dmap->get(0);
  cout<<" Finished road finding with "<< iroad.count() << " 1d roads.\n";  
  while(TMui1DRoadMapO::pointer roadptr = iroad.next())
  {
    roadptr->get()->print();
  }

  TMuiRoadMapO::iterator iroad2 = _roadmap->get(0);
  cout<<" Finished road finding with "<< iroad2.count() << " 2d roads.\n";
  while(TMuiRoadMapO::pointer roadptr = iroad2.next())
  {
    roadptr->get()->print();
  }
#endif

  _timer.get()->stop();
  return 0;
}

PHBoolean
mMuiRoadFinder1::cut_1droads()
{
  for(int iArm = 0; iArm < TMuiChannelId::kArmsTotal; iArm++)
  {
  //Now we have a list of roads tracked all the way through the
  //detector. Before we add them to the 1D Road Container
  //perform all the 1DRoad Cuts except the Duplicate cut.


    TMui1DRoadMapO::iterator iRoad1 = _road1dmap->get(iArm);

    while(TMui1DRoadMapO::pointer pRoad1 = iRoad1.next())
    {

      // Ok we found a road already, we need to know
      // whether we should save it or not
      // check road quality...
      //
#ifdef MUIOO_DEBUG
      cout<<"Checking Road: "<<endl;
      pRoad1->get()->print();
#endif
      // 1. last plane of the road
      bool LastPlaneOK = false;
      if ( pRoad1->get()->get_depth() >= _mod_par->get_min_last_gap_1d() )
      {
        LastPlaneOK = true;
      }

      // 2. positon at reference plane
      bool VertexCutOK = false;

      TMutFitPar roadfit = pRoad1->get()->get_fit_par();
     if(pRoad1->get()->get_orientation()==0)
      {
        if ( fabs(roadfit.get_y()) <= _mod_par->get_max_yref_1d() )
        {
          VertexCutOK = true;
        }
      } else {
        if ( fabs(roadfit.get_x()) <= _mod_par->get_max_xref_1d() )
        {
          VertexCutOK = true;
        }
      }

     // 3. number of gaps contain hits
     bool FiredGapsOK = false;
#ifdef MUIOO_DEBUG
     cout<<"Testing FiredGaps: G "
         << pRoad1->get()->get_gapbit()
       << " N " << pRoad1->get()->get_numfired() <<endl;
#endif
      if ( pRoad1->get()->get_numfired() > _mod_par->get_min_fired_gaps())
      {
        FiredGapsOK = true;
      }

      // Here is our criteria for a candidate road
      // at this stage.
      // i.e. 1) There are at least two gaps containing hits
      //      2) LastPlane of the road passes last plane cut
      //      3) Reference position passes vertex cut


      // At this point we are only cutting exact duplicates
      // resulting from the multiple seed loops

      if ( LastPlaneOK && VertexCutOK&& FiredGapsOK ){
        int NumHits1 = 0;
        int NumHits2 = 0;
        TMui1DRoadMapO::iterator iRoad2 = iRoad1;
        while(TMui1DRoadMapO::pointer pRoad2 = iRoad2.next()){
          if(pRoad1->get()->get_orientation() != pRoad2->get()->get_orientation()){
            continue;
          }
          TMuiClusterMapO::const_key_iterator iclust1
            =pRoad1->get()->get_associated<TMuiClusterO>();
          TMuiClusterMapO::const_key_iterator iclust2
            =pRoad2->get()->get_associated<TMuiClusterO>();
          
          NumHits1 = iclust1.count();
          NumHits2 = iclust2.count();
          
          if(NumHits1!=NumHits2){
            continue;
          }
          short SharedHits1D = 0;  
          while(TMuiClusterMapO::const_pointer clust1 = iclust1.next()){
#ifdef MUIOO_DEBUG
            clust1->get()->get_key().print();
#endif
            iclust2.reset(); 
            while(TMuiClusterMapO::const_pointer clust2 = iclust2.next()){
#ifdef MUIOO_DEBUG
              clust2->get()->get_key().print();
#endif
              if(clust1->get()->get_key()==clust2->get()->get_key()) {
                SharedHits1D++;
              }
            }
          }
#ifdef MUIOO_DEBUG          
          if(pRoad1->get()->get_panel()==pRoad1->get()->get_panel()){
            cout<<"Checked for duplicate: SharedHits("<<SharedHits1D<<") "
                <<NumHits1<<"/"<<NumHits2<<endl;
            pRoad1->get()->print();
            pRoad2->get()->print();         
          }
#endif          
          if(SharedHits1D==NumHits1){
#ifdef MUIOO_DEBUG
            cout<<"Duplicate road found with "<<SharedHits1D<<" shared hits"<<endl;
            pRoad1->get()->print();
            pRoad2->get()->print();
#endif
            // Instead of erasing now.  Lets flagit by
            // setting its depth to 0.
            // Later we'll come back and erase all.
            pRoad2->get()->set_depth(0);
          }// if SharedHits1D==NumHits1
        }//while(pRoad2
      }else{//if(VtxOk ...)
#ifdef MUIOO_DEBUG        
        cout<<"Cutting 1d road: " << "LastPlaneOK "<<LastPlaneOK
            <<" VertexCutOK " << VertexCutOK
            <<" FiredGapsOK " <<FiredGapsOK
            <<endl;
#endif
        // Instead of erasing now.  Lets flagit by
        // setting its depth to 0.
        // Later we'll come back and erase all.
        pRoad1->get()->set_depth(0);
      }
      
    }//while(pRoad1
    iRoad1 = _road1dmap->get(iArm);
    while(TMui1DRoadMapO::pointer pRoad1 = iRoad1.next())
    {
      if(pRoad1->get()->get_depth()==0)
      {
        _road1dmap->erase(pRoad1->get()->get_key());
      }
    }
  }//for(iArm..
  return 0;
}

PHBoolean
mMuiRoadFinder1::find_1droads() {
  
  // We have all of the input and output parameters now, so do
  // something useful with them!

  //
  // Now find some roads!
  //

  short Count0,iPlane0, iPlane1;
  EOrient_t orient;
 int rawcountperarm[2] = {0};

  rawcountperarm[0] = _clustermap->get(0).count();
  // JLN - the index below is now fixed 07/25/03
  rawcountperarm[1] = _clustermap->get(1).count();

  for (short iArm=0; iArm<TMuiChannelId::kArmsTotal; iArm++) {

    if(rawcountperarm[iArm]>_mod_par->get_max_occupancy_per_arm())
    {
      cout<<PHWHERE<<" Skipping analysis of arm "<<iArm<<" with too many hits "<<rawcountperarm[iArm]<<endl;
      continue;
    }
    
    for (short iorient=0; iorient<TMuiChannelId::kOrientations; iorient++){

      if(iorient==0)
        orient = kHORIZ;
      else
        orient = kVERT;

      for( short iLoop=0; iLoop<_mod_par->get_num_seed_loops(); iLoop++){

        // checkout the seed gap(s) from search order
        
        iPlane0 = _mod_par->get_search_order(iLoop,0);
        iPlane1 = _mod_par->get_search_order(iLoop,1);

        //Create a temporary list of roads for this search order.
        //list<TMui1DRoad*> roadList;

        short TotalPanel0;
        if (iPlane0 == -1){
          // one-seed-gap method, by pass iPlane0
          TotalPanel0 = 1;
        } else {
          // two-seed-gap method
          TotalPanel0 = TMuiChannelId::kPanelsPerPlane;
        }

        for (short iPanel0=0; iPanel0<TotalPanel0; iPanel0++) {

          TMuiClusterMapO::iterator iclust0;
          if (iPlane0 != -1){
            iclust0 = _clustermap->get(iArm,iPlane0,iPanel0,orient);
          }
          Count0 = (iPlane0 == -1)? (1) : iclust0.count();
          
          for (short idClust0=0; idClust0<Count0; idClust0++) {
            TMuiClusterMapO::pointer clust0ptr = iclust0.next();
            //iclust0->current() will either point to the first cluster or nothing
            
            for (short iPanel1=0; 
                 iPanel1<TMuiChannelId::kPanelsPerPlane; iPanel1++) {

              TMuiClusterMapO::iterator iclust1
                     = _clustermap->get(iArm,iPlane1,iPanel1,orient);;
       
              while(TMuiClusterMapO::pointer pC1 = iclust1.next()) {
          
                if (!(pC1->get())) {
                  cout << PHWHERE << " plane " << iPlane1 
                       << " null pointer" << endl;
                } else {

                  // Found seed cluster(s), create a new 1D road
                  // object, and attach cluster pC1 (and pC0, if 2 seed gaps)

                  TMui1DRoadMapO::pointer road1D = 
                    (_road1dmap->insert_new(iArm,iPanel1,orient)).current();

                  if (!(road1D->get())) {
                    cout << PHWHERE
                         << ":failed to create 1D road!"  << endl;
                    continue;
                  }

                  if (iPlane0 != -1 ) {
                    //If Two Gap seed set vertex weight to 0
                    PHKey::associate(road1D, clust0ptr);
                    road1D->get()->set_fitweight(0, 0.0);
                  }else{
                    //If One Seed gap method use vertex in fits.
                    road1D->get()->set_fitweight(0, 1.0);
                  }

                  PHKey::associate(road1D, pC1);
#ifdef MUIOO_DEBUG
                  cout<<"Begin tracking seed: iLoop("<< iLoop
                    << ") Panel ("<<iPanel1
                    <<")"<<endl;
                  road1D->get()->print();
                  pC1->get()->print();
#endif                  
                  track_seed(road1D,iLoop,2);

                  fit1d(road1D);                  
                } // if (!pC1)
              }//for (iClust1)
            }//for(iPanel1)
          }//for(iClust0)
        }//for(iPanel0)        
      } //short (short iLoop ...)      
    } // for (short iorient ...)
  }  // for(iArm ...)
  return True;
}

PHBoolean
mMuiRoadFinder1::find_2droads() {

  // We found all possible horizontal and vertical roads in both arms
  // and now it is time to construct 2D roads base on those 1D roads
  for (short iArm=0; iArm<TMuiChannelId::kArmsTotal; iArm++)
  {
    
    TMui1DRoadMapO::iterator iRoadH = _road1dmap->get(iArm);
    while(TMui1DRoadMapO::pointer pRoadH = iRoadH.next())
    {
      if(pRoadH->get()->get_orientation()==1) continue;
      TMui1DRoadMapO::iterator iRoadV = _road1dmap->get(iArm);
      while(TMui1DRoadMapO::pointer pRoadV = iRoadV.next())
      {
        if(pRoadV->get()->get_orientation()==0) continue;


        // What are our criteria for accepting a road
        // as valid?  Check them here.
        // 1. difference of total number clusters in vRoad and in hRoad
        // 2. difference of last gaps in vRoad and hRoad
        // 3. InterSection OK or not (if X Y clusters in each gap are
        //                            in same panel or overlap region)
        // 4. Reference position check
        // 5. Chisquare check
        // 6. Last Plane check
        // 7. Duplicate road check

        // JLN - bug fix - added fabs around the depth differenece 07/25/03
        bool LastGapDeltaOK = false;
        if (abs(pRoadH->get()->get_depth() - pRoadV->get()->get_depth()) <= _mod_par->get_max_del_last_gap() )
        {
          LastGapDeltaOK = true;
        }
        // JLN - bug fix - added fabs around the hits differenece 07/25/03
        bool TotalHitsDeltaOK = false;
        if(abs(pRoadH->get()->get_nhit() - pRoadV->get()->get_nhit()) <= _mod_par->get_max_del_total_hits() )
        {
          TotalHitsDeltaOK = true;
        }

        bool IntersectionOK = true;
        TMuiClusterMapO::const_key_iterator iClustH
          =pRoadH->get()->get_associated<TMuiClusterO>();
        while(TMuiClusterMapO::const_pointer pClustH = iClustH.next())
        {
          TMuiClusterMapO::const_key_iterator iClustV
            =pRoadV->get()->get_associated<TMuiClusterO>();
          while(TMuiClusterMapO::const_pointer pClustV = iClustV.next())
          {
            if(pClustH->get()->get_plane()!=pClustV->get()->get_plane())
            {
              continue;
            }
            
            if( !(intersectionOK(pClustH,pClustV)) ){
              IntersectionOK = false;
              break;
            }
          }
        }

        bool VertexCutOK = false;
        TMutFitPar fitH = pRoadH->get()->get_fit_par();
        TMutFitPar fitV = pRoadV->get()->get_fit_par();
        
        if ( fabs(fitH.get_y()) <= _mod_par->get_max_yref_2d() &&
             fabs(fitV.get_x()) <= _mod_par->get_max_xref_2d() ) VertexCutOK = true;

        bool ChiSquareCutOK = false;
        if ( fitV.get_chi_square() <= _mod_par->get_max_xchisq()
             &&  fitH.get_chi_square() <= _mod_par->get_max_ychisq() )
        {
          ChiSquareCutOK = true;
        }
        
        bool minLastGapOK = false;
        if (pRoadH->get()->get_depth() >= _mod_par->get_min_last_gap_2d()
            || pRoadV->get()->get_depth() >= _mod_par->get_min_last_gap_2d())
        {
          minLastGapOK = true;
        }
        if ( LastGapDeltaOK &&
             TotalHitsDeltaOK &&
             IntersectionOK &&
             VertexCutOK &&
             ChiSquareCutOK &&
             minLastGapOK
             ) {
          TMuiRoadMapO::iterator iRoad = _roadmap->insert_new(iArm);
          make2d(iRoad.current(),pRoadH,pRoadV);
        }// if(pass cuts ..
      } // for ( vRoad ..
    } // for ( hRoad ..
  }//for(iArm ...
  return 0;
}

int mMuiRoadFinder1::track_seed(TMui1DRoadMapO::pointer road1dptr,
                                int iLoop, int iSearch)
{

  double Win = _mod_par->get_mui_window();
  if(road1dptr->get()->get_associated<TMuiClusterO>().count() == 1)
    Win = 2.0 * Win;
  
  //Invalid iLoop
  if(iLoop>=_mod_par->get_num_seed_loops()) return 0;

  // End this recursion loop over search planes
  if(iSearch>=_mod_par->get_search_length(iLoop)) return 0;
#ifdef MUIOO_DEBUG
  cout<<"Tracking seed: iLoop("<<iLoop<<") iSearch("<<iSearch<<")\n";
#endif
  //First get list of clusters in panel
  // and order them by their distance from seed projection
  TMuiClusterMapO::iterator iclusts =
    _clustermap->get(road1dptr->get()->get_arm(),
                     _mod_par->get_search_order(iLoop, iSearch),
                     road1dptr->get()->get_panel(),
                     road1dptr->get()->get_orientation()
                     );



  fit1d(road1dptr);
  
  // Loop over all clusters in the panel
  // and project the road to the current panel;
  // if in the window append to local copy
  // then sort by distance
  std::vector<muioo_cluster_distance> clusterList;
  muioo_cluster_distance tmp_cd;
  double clust_dist = 0.0;
  bool clustInPanel = false;

  while(TMuiClusterMapO::pointer clustptr = iclusts.next())

  {
    clust_dist = ClusterDistance(road1dptr,clustptr);
    // JLN - this was MuidWindow, but changed to use Win (as defined above to double in case of first hit) - 07/25/03
    if(clust_dist < Win)
    {
      tmp_cd.distance = clust_dist;
      tmp_cd.pCluster = clustptr;
      clusterList.push_back(tmp_cd);
    }
  }
  
  // If there's no cluster found
  // lets look for clusters in the adjacent panels, but
  // continue to consider a candidate with no hit in the
  // current panel, ie Hardware Inefficiency
#ifdef MUIOO_DEBUG
  cout<<"Found "<<clusterList.size()<<" in the current panel "<<endl;
#endif
  if(clusterList.size()>0)
  {
    clustInPanel=true;
  }

  if(!clustInPanel)
  {
    for (short iPanel=0;
         iPanel<TMuiChannelId::kPanelsPerPlane; iPanel++) {

      if(!(
         (iPanel = ((iPanel + 1)% TMuiChannelId::kPanelsPerPlane))
         ||((iPanel + TMuiChannelId::kPanelsPerPlane - 1)% TMuiChannelId::kPanelsPerPlane)
        )) continue; //Skip non-adjacent panels
      
      TMuiClusterMapO::iterator iclustsNA =
        _clustermap->get(road1dptr->get()->get_arm(),
                         _mod_par->get_search_order(iLoop, iSearch),
                         iPanel,
                         road1dptr->get()->get_orientation()
                         );

      while(TMuiClusterMapO::pointer clustptr = iclustsNA.next())
      {

        if(!clustptr)
        {
          cout <<PHWHERE << ": null "
               << "cluster pointer"   << endl;
          break;
        }else{

          // Get the displacement of cluster pc to road1D
          clust_dist  = ClusterDistance(road1dptr,clustptr);

          // Consider all of the clusters within the window
          // JLN - again MuidWindow -> Win - 07/25/03
          if ( clust_dist < Win )
          {
            tmp_cd.pCluster=clustptr;
            tmp_cd.distance = clust_dist;
            clusterList.push_back(tmp_cd);
          }
        }//if(!pC)
      } //while(iClust)
    }//if(iPanel ...
  }//  if(!clustInPanel)
  
  //Let's sort the list by the distance from projection
  sort(clusterList.begin(), clusterList.end(), muioo_cluster_distance_sort());

  // First case: No clusters in any panel so go to next search gap
  if((!clustInPanel&& (clusterList.size()==0)))
  {
    track_seed(road1dptr , iLoop, iSearch+1);
  }else if(!clustInPanel)
  {
    //Second case: no clust in current panel but in adjacent
    // track a copy without clusters
    TMui1DRoadMapO::pointer new1droadptr = (_road1dmap->insert_new(road1dptr)).current();
#ifdef MUIOO_DEBUG
    cout<<"Copied 1DRoad\n";
    cout<<"Original:\n";
    road1dptr->get()->print();
    cout<<"Copy:\n";
    new1droadptr->get()->print();
#endif    
    track_seed(new1droadptr, iLoop, iSearch+1);
  }

  // In any case,
  // Loop over all clusters up to maxClustersPerGapSearch
  // Create a copy (if not last)
  // Add cluster to copy and track_seed
  int newRoadCount = 0;
  vector<muioo_cluster_distance>::iterator lastEntry =
    clusterList.end();
  lastEntry--;
  for(vector<muioo_cluster_distance>::iterator iClustCand = clusterList.begin();
      iClustCand != clusterList.end(); iClustCand++){
    if(iClustCand == lastEntry ||
       (newRoadCount == _mod_par->get_max_clusters_per_gap_search()-1)){
      PHKey::associate(road1dptr,(*iClustCand).pCluster);
      track_seed(road1dptr, iLoop, iSearch+1);
      break;
    }else{
      TMui1DRoadMapO::pointer new1droadptr = (_road1dmap->insert_new(road1dptr)).current();
#ifdef MUIOO_DEBUG
      cout<<"Copied 1DRoad\n";
      cout<<"Original:\n";
      road1dptr->get()->print();
      cout<<"Copy:\n";
      new1droadptr->get()->print();
      return 0;
#endif
      PHKey::associate(new1droadptr,(*iClustCand).pCluster);
      track_seed(road1dptr, iLoop, iSearch+1);
      newRoadCount++;
    }
  }
  clusterList.clear();
  return 0;
}

int mMuiRoadFinder1::make2d(TMuiRoadMapO::pointer cRoad,
                           TMui1DRoadMapO::pointer pRoadH,
                           TMui1DRoadMapO::pointer pRoadV)
{

  TMutFitPar fitH = pRoadH->get()->get_fit_par();
  TMutFitPar fitV = pRoadV->get()->get_fit_par();

  TMuiClusterMapO::key_iterator iClustH
    =pRoadH->get()->get_associated<TMuiClusterO>();
  while(TMuiClusterMapO::pointer pClust = iClustH.next())
  {
    PHKey::associate(pClust, cRoad);
  }
  TMuiClusterMapO::key_iterator iClustV
    =pRoadV->get()->get_associated<TMuiClusterO>();
  while(TMuiClusterMapO::pointer pClust = iClustV.next())
  {
    PHKey::associate(pClust,cRoad);
  }
  cRoad->get()->set_nhit(pRoadH->get()->get_nhit() + pRoadV->get()->get_nhit());
  //Combine the FitPar of the 1D roads
  TMutFitPar fit2d;
  fit2d.set_x(fitV.get_x());
  fit2d.set_y(fitH.get_y());
  fit2d.set_z(fitH.get_z());
  fit2d.set_dxdz(fitV.get_dxdz());
  fit2d.set_dydz(fitH.get_dydz());
  if(pRoadV->get()->get_depth()>pRoadH->get()->get_depth())
  {
    cRoad->get()->set_depth(pRoadV->get()->get_depth());
    fit2d.set_z_end(fitV.get_z_end());
  }else{
    cRoad->get()->set_depth(pRoadH->get()->get_depth());
    fit2d.set_z_end(fitH.get_z_end());
  }
  if(abs(fitH.get_z_begin())<abs(fitV.get_z_begin()))
  {
    fit2d.set_z_begin(fitH.get_z_begin());
  }else{
    fit2d.set_z_begin(fitV.get_z_begin());
  }
  int vDOF = pRoadV->get()->get_nhit() + (int)(pRoadV->get()->get_fitweight(0)) - 2;
  int hDOF = pRoadH->get()->get_nhit() + (int)(pRoadH->get()->get_fitweight(0)) - 2 ;
  cRoad->get()->set_freedom(vDOF+hDOF);
  // JLN - fix bug - should be OR of H and V
  cRoad->get()->set_gapbit(pRoadH->get()->get_gapbit() | (pRoadV->get()->get_gapbit() << 5));
  Float_t chisq = 0.0;
  if((vDOF + hDOF) >=0) {

    chisq = (fitV.get_chi_square()*vDOF
             +fitH.get_chi_square()*hDOF
             )/(vDOF + hDOF);
  }
  fit2d.set_chi_square(chisq);

  cRoad->get()->set_fit_par(fit2d);
  // previously not set parameters
  cRoad->get()->set_road_quality(chisq);
  UShort_t gapmatch = pRoadH->get()->get_gapbit() & 
    pRoadV->get()->get_gapbit();
  if (gapmatch != 0) {
    cRoad->get()->set_max_hit_plane(2);
  }
  else {
    cRoad->get()->set_max_hit_plane(1);
  }

  return 0;
}

int mMuiRoadFinder1::fit1d(TMui1DRoadMapO::pointer road1dptr)
{

  float z[mMuiRoadFinder1Par::NUMFITPOINTS+1]={0.0};
  float x[mMuiRoadFinder1Par::NUMFITPOINTS+1]={0.0};
  float w[mMuiRoadFinder1Par::NUMFITPOINTS+1]={0.0};
  int npoints = 0;
  float slope = 0.0,
    intercept = 0.0,
    chisq = 0.0,
    sigma_slope = 0.0,
    sigma_intercept = 0.0,
    hit_zmin = 1E10,
    hit_zmax = 0.0;

  z[npoints] = _mod_par->get_zvert();
  x[npoints] = (road1dptr->get()->get_orientation()==0)
                 ? _mod_par->get_yvert() : _mod_par->get_xvert();
  w[npoints] = _mod_par->get_weight_par_1d(0);
  npoints++;
  
  TMuiClusterMapO::const_key_iterator iclust
     =road1dptr->get()->get_associated<TMuiClusterO>();
#ifdef MUIOO_DEBUG
  cout << "mMuiRoadFinder1::fit1d" <<endl;
  road1dptr->get()->print();
  cout << "Associated cluster count("
       << iclust.count() <<")"<<endl;
#endif 
  road1dptr->get()->set_nhit(iclust.count());
  
  while(TMuiClusterMapO::const_pointer clustptr = iclust.next())
  {
    PHPoint centroid = clustptr->get()->get_centroidpos();
    PHPoint sigma = clustptr->get()->get_centroidsigma();
    road1dptr->get()->set_gapbit( road1dptr->get()->get_gapbit() | (0x1<<clustptr->get()->get_plane()));

    if(road1dptr->get()->get_depth() < clustptr->get()->get_plane())
    {
      road1dptr->get()->set_depth(clustptr->get()->get_plane());
    }
      
    z[npoints]  = centroid.getZ();
    if(abs(z[npoints])<abs(hit_zmin))
    {
      hit_zmin = z[npoints];
    }
    if(abs(z[npoints])>abs(hit_zmax))
    {
      hit_zmax = z[npoints];
    }
    if (road1dptr->get()->get_orientation() == 0) {
      x[npoints] = centroid.getY();
      w[npoints] = 1.0/(sigma.getY()*sigma.getY());
    } else {
      x[npoints] = centroid.getX();
      w[npoints] = 1.0/(sigma.getX()*sigma.getX());
    }
    npoints++;
    if(npoints>6) break;
  }

  // DS: kludge to remove vertex point from fit if we have enough points to make a line without it
  if (npoints > 2) {
    for (int ipoint = 1; ipoint<npoints; ipoint++) {
      x[ipoint-1] = x[ipoint];
      z[ipoint-1] = z[ipoint];
      w[ipoint-1] = w[ipoint];
    }
    npoints--;
  }
  // end kludge

  int status = utiLineFit(z,x,w,npoints,&slope,&intercept,&chisq,
                          &sigma_slope,&sigma_intercept);
  if(status!= 0 )
  {
    cout<<PHWHERE<<" Unable to fit MuID road." << endl;
    return -1;
  }
  double dydz, dxdz, xint, yint;
  if (road1dptr->get()->get_orientation() == 0)
  {
    dydz = slope;
    dxdz = 0.0;
    yint = intercept;
    xint = 0.0;
  } else {
    dydz = 0.0;
    dxdz = slope;
    yint = 0.0;
    xint = intercept;
  }
  TMutFitPar fitpar1d(xint,yint,0.0,dxdz,dydz,chisq);
  fitpar1d.set_z_begin(hit_zmin);
  fitpar1d.set_z_end(hit_zmax);
  road1dptr->get()->set_fit_par(fitpar1d);
  
  return 0;
}

double mMuiRoadFinder1::ClusterDistance(TMui1DRoadMapO::pointer road1dptr, TMuiClusterMapO::pointer clustptr)
{
  PHVector dirV = road1dptr->get()->get_fit_par().get_tangent();
  dirV.normalize();
  PHPoint panelX
  = TMuiGeometry::Geom()->FindIntersection(
                                           road1dptr->get()->get_arm(),
                                           clustptr->get()->get_plane(),
                                           clustptr->get()->get_panel(),
                                           road1dptr->get()->get_fit_par().get_point(),
                                           dirV
                                           );

  double distance = (road1dptr->get()->get_orientation()==0)
    ? abs(panelX.getY() - clustptr->get()->get_centroidpos().getY())
    : abs(panelX.getX() - clustptr->get()->get_centroidpos().getX());

  return distance;
}

void mMuiRoadFinder1::set_interface_ptrs(PHCompositeNode* top_node)
{
  // module runtime parameters
  try {
    _mod_par = TMutNode<mMuiRoadFinder1Par>::find_node(top_node,"mMuiRoadFinder1Par");
  }
  catch(std::exception& e){
    MUIOO::TRACE(e.what());
    return;
  }

  //Find the TMuiRoadO IOC
  try{
    _roadmap = TMutNode<TMuiRoadMapO>::find_node(top_node,"TMuiRoadMapO");
  } catch(std::exception& e){
MUIOO::TRACE(e.what());
    return;
//    TMutNode<TMuiRoadMapO>::new_node(muiNode,"TMuiRoadMapO");
//    _roadmap = TMutNode<TMuiRoadMapO>::find_node(muiNode,"TMuiRoadMapO");
//    _roadmap->make_persistent(dstNode,"TMuiRoadO");
  }

  //Find the TMui1DRoadO IOC
  try{
    _road1dmap = TMutNode<TMui1DRoadMapO>::find_node(top_node,"TMui1DRoadMapO");
  } catch(std::exception& e){
MUIOO::TRACE(e.what());
    return;
//    TMutNode<TMui1DRoadMapO>::new_node(muiNode,"TMui1DRoadMapO");
//    _road1dmap = TMutNode<TMui1DRoadMapO>::find_node(muiNode,"TMui1DRoadMapO");
//    _road1dmap->make_persistent(dstNode,"TMui1DRoadO");
  }

  //Find the new TMuiClusterO IOC
  try{
    _clustermap = TMutNode<TMuiClusterMapO>::find_node(top_node,"TMuiClusterMapO");
  } catch(std::exception& e){
MUIOO::TRACE(e.what());
    return;
 //   TMutNode<TMuiClusterMapO>::new_node(muiNode,"TMuiClusterMapO");
 //   _clustermap = TMutNode<TMuiClusterMapO>::find_node(muiNode,"TMuiClusterMapO");
 //   _clustermap->make_persistent(dstNode,"TMuiClusterO");
  }
  
}

bool mMuiRoadFinder1::intersectionOK(TMuiClusterMapO::const_pointer clustH,
                                     TMuiClusterMapO::const_pointer clustV) const
{
  TMuiGeometry* geom = TMuiGeometry::Geom();
  // Is the intersection point inside the panel(s) containing the H and
  // V clusters?
  TMuiPanelGeo* ph = geom->getPanel(clustH->get()->get_arm(),
                                    clustH->get()->get_plane(),
                                    clustH->get()->get_panel());
  TMuiPanelGeo* pv = geom->getPanel(clustV->get()->get_arm(),
                                    clustV->get()->get_plane(),
                                    clustV->get()->get_panel());

  PHPoint hpoint = clustH->get()->get_centroidpos();
  PHPoint vpoint = clustV->get()->get_centroidpos();

  PHPoint h1_local = ph->TransformToPanel(hpoint);
  PHPoint v1_local = ph->TransformToPanel(vpoint);

  PHPoint h2_local = pv->TransformToPanel(hpoint);
  PHPoint v2_local = pv->TransformToPanel(vpoint);

  return (ph->IsInPanel(v1_local.getX(), h1_local.getY(), 0.0) &&
          pv->IsInPanel(v2_local.getX(), h2_local.getY(), 0.0)   );
  
}

int mMuiRoadFinder1::flag_golden()
{
  //Loop over arms
  for (short iArm=0; iArm<TMuiChannelId::kArmsTotal; iArm++)
  {
    // Reject ghost by grouping the roads
    vector<TMuiRoadOGroup*> pGroup;
    vector<TMuiRoadOGroup*>::iterator ig;
    
    TMuiRoadOGroup* iGroup;
    TMuiRoadMapO::iterator iRoad = _roadmap->get(iArm);
    while(TMuiRoadMapO::pointer road2D = iRoad.next())
    {
      bool InGroup = false;
      iGroup=NULL;
      for(ig=pGroup.begin();ig!=pGroup.end();ig++){
        if ((*ig)->IsGroup(road2D)){
          if(!InGroup)
          {
            //This is the first group associated with the road
            (*ig)->AttachRoad(road2D);
            iGroup = (*ig);
          }else{
            // This road has already been associated with a road
            // Instead of adding it to two groups lets merge this into the first
            // then remove the first group
            // cout<<"Merging MuidGroups ..."<<endl;
            iGroup->MergeGroup(*ig);
            TMuiRoadOGroup* removedGroup = (*ig);
            delete removedGroup;  // Delete the new TMuiRoadOGroup object
            pGroup.erase(ig);  // Remove the Group pointer from the list.
            ig--; // Decrement the iterator since erase incremented
          }
          InGroup = true;
        }
      }
      if( !InGroup){
        // This road doesn't belong to any group, create a new group
        float MutrZ;
        if(iArm == 0){
          MutrZ = _mod_par->get_mut_z_south();
        } else {
          MutrZ = _mod_par->get_mut_z_north();
        }
        iGroup = new TMuiRoadOGroup(iArm,MutrZ,_mod_par->get_mut_window(),_mod_par->get_mui_window());
        iGroup->AttachRoad(road2D);
        pGroup.push_back(iGroup);
      }
    }
    
#ifdef MUIOO_DEBUG
  cout<<PHWHERE<<" Mark Golden Roads"<< endl;
#endif

    // set group index and golden mark for each road object
    // The golden flags are saved in the TMui2DRoads
    short ngroup = pGroup.size();
    for(short i=0; i<ngroup; i++){
      iGroup = pGroup[i];
      iGroup->SetGroup(i);
    }

#ifdef MUIOO_DEBUG
    cout<<PHWHERE<<" Destroy Group Objects"<< endl;
#endif

    // destroy the all group objects
    for(short i=0; i<ngroup; i++){
      delete pGroup[i];
    }

  }//for(iArm...
  return 0;
}

---