---

mMuiClusterFinder.cxx

Go to the documentation of this file.

#include <iostream>
#include "PHNode.h"
#include "PHIODataNode.h"
#include "PHNodeIterator.h"

#include "mMuiClusterFinder.h"

typedef PHIODataNode<PHTable> TableNode_t;
// STL hit list iterator
typedef list<TMuiHitO>::iterator hit_list_iterator;

mMuiClusterFinder::mMuiClusterFinder():
  _timer( PHTimeServer::get()->insert_new("mMuiClusterFinder") )
  {}

mMuiClusterFinder::~mMuiClusterFinder(){}

PHBoolean
mMuiClusterFinder::event(PHCompositeNode *top_node)
{
  _timer.get()->restart();
  try {     
    // Reset IOC pointers
    //
    set_interface_ptrs(top_node);
    
    // Associate groups of contiguous hits with 
    // TMutClus objects
    //
    find_clusters();
        
  } catch(std::exception& e) {
    MUIOO::TRACE(e.what());
    return False;
  }
  _timer.get()->stop();  
  return True;  
}

void
mMuiClusterFinder::set_interface_ptrs(PHCompositeNode* top_node)
{
   // module runtime parameters
  _mod_par = TMutNode<mMuiClusterFinderPar>::find_node(top_node,"mMuiClusterFinderPar");

  //Find the new TMuiHitO objects
  _muihitmap = TMutNode<TMuiHitMapO>::find_node(top_node,"TMuiHitMapO");

  //Find the new TMuiClusterO objects
  _muiclustermap = TMutNode<TMuiClusterMapO>::find_node(top_node,"TMuiClusterMapO");

  return;
}

/* Find clusters from raw hits */
void
mMuiClusterFinder::find_clusters()
{
  // ===================================================================
  // Find clusters and fill the cluster container.
  // ===================================================================
  
  // For now just one hit per cluster
  // Let's make it optional though, to merge isolated neighboring hits into single clusters
  for (short arm=0; arm<TMuiChannelId::kArmsTotal; arm++) {    
    for (short plane=0; plane<TMuiChannelId::kPlanesPerArm; plane++) {
      for (short panel=0; panel<TMuiChannelId::kPanelsPerPlane; panel++) {

        //Lets get a pointer to this panel's geometry object
        TMuiChannelId pPanelId(arm,plane,panel);
        TMuiPanelGeo* fpMuiPanelGeo   = TMuiGeometry::Geom()->getPanel(pPanelId);
        
        for (short orient=0; orient<TMuiChannelId::kOrientations; orient++) {
          
          TMuiHitMapO::iterator hit_iter = _muihitmap->get(arm,plane,panel,orient);
          while(TMuiHitMapO::pointer hit_ptr = hit_iter.next()) {
            if (_mod_par->get_clustering_mode() == mMuiClusterFinderPar::NONE) {
              // no clustering, just make each hit into its own cluster
              make_new_cluster(hit_ptr);
            }
            else {
              // iterator to all clusters with same location as hit
              TMuiClusterMapO::iterator clus_iter = _muiclustermap->get(arm, plane, panel, orient);
              bool hit_is_added = false;
              while(TMuiClusterMapO::pointer clus_ptr = clus_iter.next()){
                if (!hit_is_added || _mod_par->get_clustering_mode() == mMuiClusterFinderPar::NORMAL)
                  { // in the 'NORMAL' mode, we can try to add the hit to as many clusters as possible
                    // attempt to add hit to cluster
                    hit_is_added = associate_hit(clus_ptr,hit_ptr); 
                  }
              }
              // hit was not associated with existing cluster so make
              // a new cluster and associate this hit
              if(!hit_is_added){
                make_new_cluster(hit_ptr);
              }
            }            
          }//loop over hits in panel/orient

          // Now, the clustering is done. Let's calculate the positions if we found any
          TMuiClusterMapO::iterator clus_iter = _muiclustermap->get(arm, plane, panel, orient);
          while(TMuiClusterMapO::pointer clus_ptr = clus_iter.next()){
            calc_position(clus_ptr, fpMuiPanelGeo);
          }

        }//orient
      }//panel
    }//plane
  }//arm
  
  return;
}

bool
mMuiClusterFinder::associate_hit(TMuiClusterMapO::pointer clus_ptr, 
                                 TMuiHitMapO::pointer hit_ptr)
{  
  // iterator to hits associated with this cluster
  //
  TMuiHitMapO::key_iterator hit_iter = 
    clus_ptr->get()->get_associated<TMuiHitO>();
  
  // loop over associated hits
  //
  while(TMuiHitMapO::pointer clus_hit_ptr = hit_iter.next()){

    UShort_t distance = std::abs(clus_hit_ptr->get()->get_twopack() - 
                                 hit_ptr->get()->get_twopack()); 
    if(distance == 1) {
      // strips are adjacent so associate hit with cluster
      //      
      PHKey::associate(hit_ptr,clus_ptr);
      return true;
    }
  }
  // hit was not associated return false
  //
  return false;
}

void
mMuiClusterFinder::calc_position(TMuiClusterMapO::pointer clus_ptr,
                                 TMuiPanelGeo *fpMuiPanelGeo)
{  
  // iterator to hits associated with this cluster
  //
  TMuiHitMapO::key_iterator hit_iter = 
    clus_ptr->get()->get_associated<TMuiHitO>();
  
  // loop over associated hits and get their positions
  //
  PHVector clpos;
  PHVector clerr;
  int clnhit = 0;
  EOrient_t orient = kHORIZ;
  while(TMuiHitMapO::pointer clus_hit_ptr = hit_iter.next()){
    
    TMuiChannelId pTPID(clus_hit_ptr->get()->get_arm(),
                        clus_hit_ptr->get()->get_plane(),
                        clus_hit_ptr->get()->get_panel(),
                        (EOrient_t)clus_hit_ptr->get()->get_orientation(),
                        clus_hit_ptr->get()->get_twopack());

    orient = (EOrient_t)clus_hit_ptr->get()->get_orientation();
    TMuiTwoPackGeo* fpMuiTwoPackGeo = TMuiGeometry::Geom()->getTwoPack(pTPID);
    
    //Now we need to find the centroid of this tube
    if (fpMuiPanelGeo && fpMuiTwoPackGeo) {
      float x, y, z, sx, sy, sz;
      
      fpMuiTwoPackGeo->CenterPos(x, y, z);
      fpMuiTwoPackGeo->CenterSigma(sx, sy, sz);
      PHVector pos(x, y, z);
      PHVector sigma(sx, sy, sz);
      PHVector fCentroidPos = fpMuiPanelGeo->TransformToGlobal(pos);
      PHVector fCentroidErr = fpMuiPanelGeo->RotateToGlobal(sigma);

      clpos = clpos + fCentroidPos;
      clerr = clerr + fCentroidErr;
      clnhit++;
    }
  }

  if (clnhit > 0)
    {
      clpos = clpos * (1.0/clnhit);
      // we should scale the error also
      // if vert., the values/errors are the same for the two-packs in y,z 
      // if horiz., the values/errors are the same for the two-packs in x,z 
      //
      if (orient == kHORIZ)
        { // horizontal, don't scale y-value
          clerr.setX(clerr.getX() *( 1.0/clnhit) );
          clerr.setZ(clerr.getZ() *( 1.0/clnhit) );
        }
      else
        { // vertical, don't scale x-value
          clerr.setY(clerr.getY() *( 1.0/clnhit) );
          clerr.setZ(clerr.getZ() *( 1.0/clnhit) );
        }
      clus_ptr->get()->set_centroidpos(PHPoint(clpos));
      clus_ptr->get()->set_centroidsigma(PHPoint(clerr));
      clus_ptr->get()->set_size(clnhit);

      // figure out what the start and end points are also
      // using the centroid and sigma info
      //
      PHPoint startpoint = PHPoint(clpos);
      PHPoint endpoint = startpoint;
      // from sigma to total range of uncertainty: sqrt(12) conv. factor
      // and then we want half the length
      double convfactor = sqrt(12.0)*0.5;
      PHPoint diff = PHPoint(clerr)*convfactor;

      //
      if (orient == kHORIZ)
        { // horizontal
          startpoint.setX(startpoint.getX() - diff.getX());
          endpoint.setX(endpoint.getX() + diff.getX());
        }
      else
        { // vertical
          startpoint.setY(startpoint.getY() - diff.getY());
          endpoint.setY(endpoint.getY() + diff.getY());
        }
      // and store the values
      clus_ptr->get()->set_coord_end(endpoint);
      clus_ptr->get()->set_coord_begin(startpoint);
    }
  else
    { // couldn't locate coordinates of any hits.
      // That must be bad. We might as well remove this cluster
      _muiclustermap->erase(clus_ptr->get()->get_key());
    }
  return;
}

void
mMuiClusterFinder::make_new_cluster(TMuiHitMapO::pointer hit_ptr)  
{
  // make a new cluster with same location as hit
  //
  TMuiClusterMapO::iterator clus_iter = 
    _muiclustermap->insert_new(hit_ptr->get()->get_arm(),
                               hit_ptr->get()->get_plane(),
                               hit_ptr->get()->get_panel(),
                               hit_ptr->get()->get_orientation());
  
  // Associate new cluster with hit
  //
  PHKey::associate(hit_ptr,clus_iter.current());
}

void mMuiClusterFinder::apply_cluster_cuts()
{
  TMuiClusterMapO::iterator clus_iter = _muiclustermap->range();
  while(TMuiClusterMapO::pointer clus_ptr = clus_iter.next()){
    
    // If cluster contains less than the required number of
    // strips remove from map
    //
    if(clus_ptr->get()->get_size() < _mod_par->get_min_cluster_width() ) {
      _muiclustermap->erase(clus_ptr->get()->get_key());
      continue;
    }
    
    if(clus_ptr->get()->get_size() > _mod_par->get_max_cluster_width() ) {
      _muiclustermap->erase(clus_ptr->get()->get_key());
      continue;
    }
  }
}

---