PHTpcSeedFinder.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTpcSeedFinder.cc

#include "PHTpcSeedFinder.h"
#include "PHTpcTrackerUtil.h"

#include "externals/kdfinder.hpp"  // for TrackCandidate, Circle, find_tr...

#include <phool/PHLog.h>

#include <log4cpp/CategoryStream.hh>  // for CategoryStream

class TrkrClusterContainer;
class TrkrHitSetContainer;

float round(float var)
{
  return (int(var * 100.0) / 100.0);
}

PHTpcSeedFinder::PHTpcSeedFinder()
  : mMaxDistance1(3.0)
  , mTripletAngle1(M_PI / 8)
  , mMinHits1(10)
  , mMaxDistance2(6.0)
  , mTripletAngle2(M_PI / 8)
  , mMinHits2(5)
  , mNThreads(1)
  , mRemoveLoopers(false)
  , mMinLooperRadius(10.0)
  , mMaxLooperRadius(70.0)
{
}

std::vector<kdfinder::TrackCandidate<double>*> PHTpcSeedFinder::findSeeds(TrkrClusterContainer* cluster_map, TrkrHitSetContainer* hitsets, double B)
{
  LOG_WARN_IF("tracking.PHTpcSeedFnder.findSeeds", !cluster_map) << __FILE__ << "," << __LINE__ << " Can't find node TRKR_CLUSTER";
  std::vector<kdfinder::TrackCandidate<double>*> result;
  if (!cluster_map)
  {
    return result;
  }
  LOG_WARN_IF("tracking.PHTpcSeedFnder.findSeeds", !hitsets) << __FILE__ << "," << __LINE__ << " Can't find node TRKR_HITSET";
  if (!hitsets)
  {
    return result;
  }

  //----- convert clusters to kdhits -----
  std::vector<std::vector<double> > kdhits(PHTpcTrackerUtil::convert_clusters_to_hits(cluster_map, hitsets));
  LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "imported " << kdhits.size() << " clusters";

  //----- find unfitted kdtracks -----
  std::vector<std::vector<double> > unused_hits;
  std::vector<std::vector<std::vector<double> > > kdtracks;

  kdtracks = kdfinder::find_tracks_iterative<double>(kdhits, unused_hits,
                                                     mMaxDistance1 /* max distance in cm*/, mTripletAngle1 /* triplet angle */, mMinHits1 /* min hits to keep track */,  // first iteration
                                                     mMaxDistance2, mTripletAngle2, mMinHits2,                                                                           // second iteration params
                                                     mNThreads /* nthreads */,
                                                     false);  // print stats

  LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "kdtracks: " << kdtracks.size();

  //----- create fitted track candidates -----
  result = kdfinder::get_track_candidates<double>(kdtracks, B);
  LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "initial kdtrack candidates: " << result.size();

  //----- filter out junk = bad tracks, <25 MeV Pt tracks etc.. -----
  for (auto it = result.begin(); it != result.end();)
  {
    if (!(*it)->isFitted() || (*it)->Pt() < 0.025 || (*it)->Pt() > 200.0)
    {
      (*it)->deleteHits();
      delete (*it);
      it = result.erase(it);
    }
    else
    {
      ++it;
    }
  }
  LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "kdtrack candidates after filtering: " << result.size();

  //----- filter out loopers, if enabled -----
  if (mRemoveLoopers)
  {
    for (auto it = result.begin(); it != result.end();)
    {
      kdfinder::Circle<double>* c = (*it)->getCircleFit();
      double x = c->a, y = c->b, r = c->r, xyr = std::sqrt(x * x + y * y);
      if (((xyr - r) > mMinLooperRadius) && ((xyr + r) < mMaxLooperRadius))
      {
        (*it)->deleteHits();
        delete (*it);
        it = result.erase(it);
      }
      else
      {
        ++it;
      }
    }
  }

  /*
  //----- primitive yet very fast vertex seed finding -----
  // FIXME: use it instead of Rave?
  std::vector< std::pair< std::vector<double>, std::vector<size_t> > > vertices =
      kdfinder::find_vertex_seeds<double>( result, 0, 0, 0.5, 2.0, 3, false );
    LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "track candidate vertex seeds: " << vertices.size();
  for( auto vertex: vertices ) {
    std::vector<double> pos = vertex.first;
    std::vector<size_t> ids = vertex.second;
    LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "vertex tracks: " << ids.size() << ", pos: " << pos[0] << ", " << pos[1] << ", " << pos[2];
  }
  */

  return result;
}