AxialHitQuadTreeProcessor Class Reference

A QuadTreeProcessor for TrackHits. More...

#include <AxialHitQuadTreeProcessor.h>

Inheritance diagram for AxialHitQuadTreeProcessor:

Collaboration diagram for AxialHitQuadTreeProcessor:

Classes
struct	NodeCache
	Geometry of a node that is needed in the containment check. More...
struct	ThetaSpanCache
	Quantities entering the node containment check that only depend on the hit and on the theta span of a node - shared by all nodes that have the same theta span. More...

Public Types
using	Item
	The QuadTree will only see items of this type.
using	QuadTree
	The used QuadTree.
using	XSpan
	This pair describes the span in X for a node.
using	YSpan
	This pair describes the span in Y for a node.
using	XYSpans
	This pair of spans describes the span of a node.
using	QuadTreeChildren
	Alias for the QuadTree Children.
using	CandidateReceiver
	This lambda function can be used for postprocessing.

Public Member Functions
	AxialHitQuadTreeProcessor (int lastLevel, int seedLevel, const XYSpans &ranges, PrecisionUtil::PrecisionFunction precisionFunction)
	Constructor.
	AxialHitQuadTreeProcessor (const TrackingUtilities::Vector2D &localOrigin, const YSpan &curvSpan, const TrackingUtilities::LookupTable< TrackingUtilities::Vector2D > *cosSinLookupTable)
	Constructor for the quad tree processor used in the off-origin extension.
void	drawHits (std::vector< const TrackingUtilities::CDCWireHit * > hits, unsigned int color=46) const
	Draw QuadTree node.
void	drawNode (QuadTree *node) const
	Draw QuadTree node.
void	clear ()
	Delete all the QuadTreeItems in the tree and clear the tree.
void	seed (const std::vector< const TrackingUtilities::CDCWireHit * > &datas)
	Fill in the items in the given vector.
std::vector< const TrackingUtilities::CDCWireHit * >	getAssignedItems ()
	Get items that have been assigned to the seed level The returned elements are unique even if items are assigned multiple times.
void	fill (const CandidateReceiver &candidateReceiver, int nHitsThreshold)
	Start filling the already created tree.
void	fill (const CandidateReceiver &candidateReceiver, int nHitsThreshold, float yLimit)
	Fill vector of QuadTree instances with hits.
virtual void	afterFillDebugHook (QuadTreeChildren &children)
	Override that function if you want to receive debug output whenever the children of a node are filled the first time Maybe you want to make some nice plots or statistics.
const std::map< std::pair< long, float >, std::vector< Item * > > &	getDebugInformation () const
	Return the debug information if collected.

Static Public Member Functions
static const TrackingUtilities::LookupTable< TrackingUtilities::Vector2D > &	getCosSinLookupTable ()
	Get the standard lookup table containing equally spaces unit vectors (cos, sin)
static std::vector< float >	createCurvBound (YSpan curvSpan, int lastLevel)
	Constructs an array with the curvature bounds as generated by the default bin divisions.

Protected Member Functions
bool	isLeaf (QuadTree *node) const final
	lastLevel depends on curvature of the track candidate
XYSpans	createChild (QuadTree *node, int i, int j) const final
	Return the new ranges.
bool	isInNode (QuadTree *node, const TrackingUtilities::CDCWireHit *wireHit) const final
	Check whether hit belongs to the quadtree node:
void	insertItemsInNodes (const std::vector< QuadTree * > &nodes, const std::vector< Item * > &items) final
	Insert the hits into the given nodes sharing the parts of the containment check that do not depend on the individual node.
bool	checkDerivative (QuadTree *node, const TrackingUtilities::CDCWireHit *wireHit) const
	Check derivative of the Legendre curve.
bool	checkExtremum (QuadTree *node, const TrackingUtilities::CDCWireHit *wireHit) const
	Checks whether extreme point is located within QuadTree node's ranges.
virtual XYSpans	createChild (QuadTree *node, int iX, int iY) const
	Implement that function if you want to provide a new processor.
virtual bool	isInNode (QuadTree *node, const TrackingUtilities::CDCWireHit *item) const=0
	Implement that function if you want to provide a new processor.
virtual void	insertItemsInNodes (const std::vector< QuadTree * > &nodes, const std::vector< Item * > &items)
	Insert each of the given items into every one of the given nodes it belongs to.
virtual bool	isLeaf (QuadTree *node) const
	Function which checks if given node is leaf Implemented as virtual to keep possibility of changing lastLevel values depending on region is phase-space (i.e.
int	getLastLevel () const
	Return the parameter last level.

Protected Attributes
std::unique_ptr< QuadTree >	m_quadTree
	The quad tree we work with.
std::deque< Item >	m_items
	Storage space for the items that are referenced by the quad tree nodes.
std::vector< Item * >	m_itemPtrs
	Reusable buffer with pointers to all items - used to seed the tree.
std::vector< QuadTree * >	m_childPtrs
	Reusable buffer with pointers to the children of the node currently being filled.
std::vector< QuadTree * >	m_seededTrees
	Vector of QuadTrees QuadTree instances (which are filled in the vector) cover the whole Legendre phase-space; each instance is processes independently.

Private Member Functions
void	fillGivenTree (QuadTree *node, const CandidateReceiver &candidateReceiver, int nItemsThreshold, float yLimit)
	Internal function to do the real quad tree search: fill the nodes, check which of the n*m bins we need to process further and go one level deeper.
void	createChildren (QuadTree *node, QuadTreeChildren &m_children) const
	Creates the sub node of a given node.
void	fillChildren (QuadTree *node, const std::vector< Item * > &items)
	This function is called by fillGivenTree and fills the items into the corresponding children.
void	callResultFunction (QuadTree *node, const CandidateReceiver &candidateReceiver) const
	When a node is accepted as a result, we extract a vector with the items (back transformed to AData*) and pass it together with the result node to the candidate receiver function.

Private Attributes
std::vector< ThetaSpanCache >	m_thetaSpanCaches
	Reusable buffer with the per theta span quantities - one entry per distinct theta span.
std::vector< NodeCache >	m_nodeCaches
	Reusable buffer with the per node quantities - one entry per node.
PrecisionUtil::PrecisionFunction	m_precisionFunction
	Lambda which holds resolution function for the quadtree.
TrackingUtilities::Vector2D	m_localOrigin
	Local origin on which the phase space coordinates are centered.
const TrackingUtilities::LookupTable< TrackingUtilities::Vector2D > *	m_cosSinLookupTable
	Pinned lookup table for precomputed cosine and sine values.
const double	c_curlCurv = 0.02
	The curvature above which the trajectory is considered a curler.
bool	m_twoSidedPhaseSpace
	Indicator whether the two sided phases space insertion check should be used This option should automatically split back to back tracks in the low curvature regions.
int	m_lastLevel
	The last level to be filled.
int	m_seedLevel
	The first level to be filled, effectively skip forward to this higher granularity level.
bool	m_debugOutput
	A flag to control the creation of the debug output.
std::map< std::pair< long, float >, std::vector< Item * > >	m_debugOutputMap
	The calculated debug map.

Static Private Attributes
static const int	c_maxNChildren
	Upper bound on the number of children of a node - they are kept on the stack.

Detailed Description

A QuadTreeProcessor for TrackHits.

Definition at line 31 of file AxialHitQuadTreeProcessor.h.

Member Typedef Documentation

CandidateReceiver

using CandidateReceiver

inherited

This lambda function can be used for postprocessing.

Definition at line 69 of file QuadTreeProcessor.h.

Item

using Item

inherited

The QuadTree will only see items of this type.

Definition at line 51 of file QuadTreeProcessor.h.

QuadTree

using QuadTree

inherited

The used QuadTree.

Definition at line 54 of file QuadTreeProcessor.h.

QuadTreeChildren

using QuadTreeChildren

inherited

Alias for the QuadTree Children.

Definition at line 66 of file QuadTreeProcessor.h.

XSpan

using XSpan

inherited

This pair describes the span in X for a node.

Definition at line 57 of file QuadTreeProcessor.h.

XYSpans

using XYSpans

inherited

This pair of spans describes the span of a node.

Definition at line 63 of file QuadTreeProcessor.h.

YSpan

using YSpan

inherited

This pair describes the span in Y for a node.

Definition at line 60 of file QuadTreeProcessor.h.

Constructor & Destructor Documentation

AxialHitQuadTreeProcessor() [1/2]

AxialHitQuadTreeProcessor	(	int	lastLevel,
		int	seedLevel,
		const XYSpans &	ranges,
		PrecisionUtil::PrecisionFunction	precisionFunction )

Constructor.

Definition at line 97 of file AxialHitQuadTreeProcessor.cc.

  : QuadTreeProcessor(lastLevel, seedLevel, ranges)
  , m_precisionFunction(precisionFunction)
  , m_localOrigin(0.0, 0.0)
  , m_cosSinLookupTable(&getCosSinLookupTable())
{
  m_twoSidedPhaseSpace = m_quadTree->getYMin() * m_quadTree->getYMax() < 0;
}

AxialHitQuadTreeProcessor() [2/2]

AxialHitQuadTreeProcessor	(	const TrackingUtilities::Vector2D &	localOrigin,
		const YSpan &	curvSpan,
		const TrackingUtilities::LookupTable< TrackingUtilities::Vector2D > *	cosSinLookupTable )

Constructor for the quad tree processor used in the off-origin extension.

Currently only used in zero level mode to collect hits that are in a phase space part with respect to the given point.

Definition at line 109 of file AxialHitQuadTreeProcessor.cc.

  : QuadTreeProcessor(0, 0, { {0, cosSinLookupTable->getNPoints() - 1}, curvSpan})
, m_localOrigin(localOrigin)
, m_cosSinLookupTable(cosSinLookupTable)
{
  // Never use two sided mode in off origin extension
  m_twoSidedPhaseSpace = false;
}

Member Function Documentation

afterFillDebugHook()

virtual void afterFillDebugHook ( QuadTreeChildren & children )

inlinevirtualinherited

Override that function if you want to receive debug output whenever the children of a node are filled the first time Maybe you want to make some nice plots or statistics.

Definition at line 392 of file QuadTreeProcessor.h.

      {
        if (not m_debugOutput) return;
        for (QuadTree& childNode : children) {
          if (childNode.getLevel() != getLastLevel()) continue; // Only write the lowest level
          //m_debugOutputMap[ {childNode.getXMean(), childNode.getYMean()}] = childNode.getItems();
        }
      }

callResultFunction()

void callResultFunction ( QuadTree * node, const CandidateReceiver & candidateReceiver ) const

inlineprivateinherited

When a node is accepted as a result, we extract a vector with the items (back transformed to AData*) and pass it together with the result node to the candidate receiver function.

Definition at line 301 of file QuadTreeProcessor.h.

      {
        const std::vector<Item*>& foundItems = node->getItems();
        std::vector<AData*> candidate;
        candidate.reserve(foundItems.size());
 
        for (Item* item : foundItems) {
          item->setUsedFlag();
          candidate.push_back(item->getPointer());
        }
 
        candidateReceiver(candidate, node);
      }

checkDerivative()

bool checkDerivative ( QuadTree * node, const TrackingUtilities::CDCWireHit * wireHit ) const

protected

Check derivative of the Legendre curve.

Parameters

node	QuadTree node
wireHit	pointer to the hit to check

Returns

returns true in cases:

- positive derivative and no extremum in the node's ranges or

- extremum located in the node's ranges

returns false in other cases (namely negative derivative

Definition at line 396 of file AxialHitQuadTreeProcessor.cc.

{
  const Vector2D& pos2D = wireHit->getRefPos2D() - m_localOrigin;
 
  long thetaMin = node->getXMin();
  long thetaMax = node->getXMax();
 
  const Vector2D& thetaVecMin = m_cosSinLookupTable->at(thetaMin);
  const Vector2D& thetaVecMax = m_cosSinLookupTable->at(thetaMax);
 
  float rMinD = thetaVecMin.cross(pos2D);
  float rMaxD = thetaVecMax.cross(pos2D);
 
  // Does not really make sense...
  if ((rMinD > 0) && (rMaxD * rMinD >= 0)) return true;
  if ((rMaxD * rMinD < 0)) return true;
  return false;
}

checkExtremum()

bool checkExtremum ( QuadTree * node, const TrackingUtilities::CDCWireHit * wireHit ) const

protected

Checks whether extreme point is located within QuadTree node's ranges.

Parameters

node	QuadTree node
wireHit	hit to check

Returns

true or false

Definition at line 415 of file AxialHitQuadTreeProcessor.cc.

{
  const double& driftLength = wireHit->getRefDriftLength();
  const Vector2D& pos2D = wireHit->getRefPos2D() - m_localOrigin;
  double r2 = pos2D.normSquared() - driftLength * driftLength;
 
  // get left and right borders of the node
  long thetaMin = node->getXMin();
  long thetaMax = node->getXMax();
 
  const Vector2D& thetaVecMin = m_cosSinLookupTable->at(thetaMin);
  const Vector2D& thetaVecMax = m_cosSinLookupTable->at(thetaMax);
 
  if (not pos2D.isBetween(thetaVecMin, thetaVecMax)) return false;
 
  // compute Legendre curves at the position
  double r = pos2D.norm();
  float rRight = r - driftLength;
  float rLeft = r + driftLength;
 
  // get top and bottom borders of the node
  float rMin = node->getYMin() * r2 / 2;
  float rMax = node->getYMax() * r2 / 2;
 
  bool crossesRight = (rMin - rRight) * (rMax - rRight) < 0;
  bool crossesLeft = (rMin - rLeft) * (rMax - rLeft) < 0;
  return crossesRight or crossesLeft;
}

clear()

void clear ( )

inlineinherited

Delete all the QuadTreeItems in the tree and clear the tree.

Definition at line 102 of file QuadTreeProcessor.h.

      {
        m_seededTrees.clear();
        m_quadTree->clearChildren();
        m_quadTree->clearItems();
        m_items.clear();
      }

createChild() [1/2]

AxialHitQuadTreeProcessor::XYSpans createChild ( QuadTree * node, int i, int j ) const

finalprotected

Return the new ranges.

We do not use the standard ranges for the lower levels.

Parameters

node	quadtree node
i	theta index of the child
j	rho index of the child

Returns

returns ranges of the (i;j) child

Definition at line 136 of file AxialHitQuadTreeProcessor.cc.

{
  const int nodeLevel = node->getLevel();
  const int lastLevel = getLastLevel();
  const float meanCurv = std::fabs(node->getYMax() + node->getYMin()) / 2;
 
  // Expand bins for all nodes 7 levels before the last level (for lastLevel = 12 starting at 6)
  // but only in a curvature region higher than 0.005. Lower than that use always standard.
  bool standardBinning = (nodeLevel <= lastLevel - 7) or (meanCurv <= 0.005);
 
  if (standardBinning) {
    float r1 = node->getYLowerBound(j);
    float r2 = node->getYUpperBound(j);
    long theta1 = node->getXLowerBound(i);
    long theta2 = node->getXUpperBound(i);
 
    // Standard bin division
    return XYSpans({theta1, theta2}, {r1, r2});
  }
 
  // Non-standard binning
  // For level 6 to 7 only expand 1 / 4, for higher levels expand  1 / 8.
  // (assuming last level == 12)
  if (nodeLevel < lastLevel - 5) {
    float r1 = node->getYLowerBound(j) - node->getYBinWidth(j) / 4.;
    float r2 = node->getYUpperBound(j) + node->getYBinWidth(j) / 4.;
 
    // long extension = pow(2, lastLevel - nodeLevel) / 4; is same as:
    long extension = pow(2, lastLevel - nodeLevel - 2);
 
    long theta1 = node->getXLowerBound(i) - extension;
    if (theta1 < 0) theta1 = 0;
 
    long theta2 = node->getXUpperBound(i) + extension;
    if (theta2 >= m_cosSinLookupTable->getNPoints()) {
      theta2 = m_cosSinLookupTable->getNPoints() - 1;
    }
 
    return XYSpans({theta1, theta2}, {r1, r2});
  } else {
    float r1 = node->getYLowerBound(j) - node->getYBinWidth(j) / 8.;
    float r2 = node->getYUpperBound(j) + node->getYBinWidth(j) / 8.;
 
    // long extension = pow(2, lastLevel - nodeLevel) / 8; is same as
    long extension = pow(2, lastLevel - nodeLevel - 3);
 
    long theta1 = node->getXLowerBound(i) - extension;
    if (theta1 < 0) theta1 = 0;
 
    long theta2 = node->getXUpperBound(i) + extension;
    if (theta2 >= m_cosSinLookupTable->getNPoints()) {
      theta2 = m_cosSinLookupTable->getNPoints() - 1;
    }
 
    return XYSpans({theta1, theta2}, {r1, r2});
  }
}

createChild() [2/2]

virtual XYSpans createChild ( QuadTree * node, int iX, int iY ) const

inlineprotectedvirtualinherited

Implement that function if you want to provide a new processor.

It decides which node-spans the n * m children of the node should have. It is called when creating the nodes. The two indices iX and iY tell you where the new node will be created (as node.children[iX][iY]). You can check some information on the level or the x- or y-values by using the methods implemented for node.

Returns

a XYSpan pair of a x- and a y-span that the new child should have. If you don nt want to provide custom spans, just return XYSpans(XSpan(node->getXBinBound(iX), node->getXBinBound(iX + 1)), YSpan(node->getYBinBound(iY), node->getYBinBound(iY + 1)));

Definition at line 324 of file QuadTreeProcessor.h.

      {
        AX xMin = node->getXLowerBound(iX);
        AX xMax = node->getXUpperBound(iX);
        AY yMin = node->getYLowerBound(iY);
        AY yMax = node->getYUpperBound(iY);
        return XYSpans({xMin, xMax}, {yMin, yMax});
      }

createChildren()

void createChildren ( QuadTree * node, QuadTreeChildren & m_children ) const

inlineprivateinherited

Creates the sub node of a given node.

This function is called by fillGivenTree. To calculate the spans of the children nodes the user-defined function createChiildWithParent is used.

Definition at line 264 of file QuadTreeProcessor.h.

      {
        m_children.reserve(node->getXNbins() * node->getYNbins());
        for (int i = 0; i < node->getXNbins(); ++i) {
          for (int j = 0; j < node->getYNbins(); ++j) {
            const XYSpans& xySpans = createChild(node, i, j);
            const XSpan& xSpan = xySpans.first;
            const YSpan& ySpan = xySpans.second;
            m_children.push_back(QuadTree(xSpan, ySpan, node->getLevel() + 1, node));
          }
        }
      }

createCurvBound()

std::vector< float > createCurvBound ( YSpan curvSpan, int lastLevel )

static

Constructs an array with the curvature bounds as generated by the default bin divisions.

Definition at line 65 of file AxialHitQuadTreeProcessor.cc.

{
  std::vector<YSpan> spans{{curvSpan}};
 
  std::vector<YSpan> nextSpans;
  for (int level = 1; level <= lastLevel; ++level) {
    nextSpans.clear();
    for (const YSpan& span : spans) {
      nextSpans.push_back(splitCurvSpan(span, level, lastLevel, 0));
      nextSpans.push_back(splitCurvSpan(span, level, lastLevel, 1));
    }
    spans.swap(nextSpans);
  }
 
  std::vector<float> bounds;
  for (const YSpan& span : spans) {
    bounds.push_back(span[0]);
    bounds.push_back(span[1]);
  }
 
  assert(bounds.size() == std::pow(2, lastLevel));
  return bounds;
}

drawHits()

void drawHits	(	std::vector< const TrackingUtilities::CDCWireHit * >	hits,
		unsigned int	color = 46 ) const

Draw QuadTree node.

Definition at line 444 of file AxialHitQuadTreeProcessor.cc.

{
  static int nevent(0);
 
  TCanvas* canv = new TCanvas("canv", "legendre transform", 0, 0, 1200, 600);
  canv->cd(1);
  TGraph* dummyGraph = new TGraph();
  dummyGraph->SetPoint(1, -M_PI, 0);
  dummyGraph->SetPoint(2, M_PI, 0);
  dummyGraph->Draw("AP");
  dummyGraph->GetXaxis()->SetTitle("#theta");
  dummyGraph->GetYaxis()->SetTitle("#rho");
  dummyGraph->GetXaxis()->SetRangeUser(-M_PI, M_PI);
  dummyGraph->GetYaxis()->SetRangeUser(-0.02, 0.15);
 
  for (const CDCWireHit* wireHit : hits) {
    const double& driftLength = wireHit->getRefDriftLength();
    const Vector2D& pos2D = wireHit->getRefPos2D() - m_localOrigin;
    double x = pos2D.x();
    double y = pos2D.y();
    double r2 = pos2D.normSquared() - driftLength * driftLength;
 
    TF1* concaveHitLegendre = new TF1("concaveHitLegendre", "2*([0]/[3])*cos(x) + 2*([1]/[3])*sin(x) + 2*([2]/[3])", -M_PI, M_PI);
    TF1* convexHitLegendre = new TF1("convexHitLegendre", "2*([0]/[3])*cos(x) + 2*([1]/[3])*sin(x) - 2*([2]/[3])", -M_PI, M_PI);
    concaveHitLegendre->SetLineWidth(1);
    convexHitLegendre->SetLineWidth(1);
    concaveHitLegendre->SetLineColor(color);
    convexHitLegendre->SetLineColor(color);
 
    concaveHitLegendre->SetParameters(x, y, driftLength, r2);
    convexHitLegendre->SetParameters(x, y, driftLength, r2);
    concaveHitLegendre->Draw("CSAME");
    convexHitLegendre->Draw("CSAME");
  }
//   canv->Print(Form("legendreHits_%i.root", nevent));
//   canv->Print(Form("legendreHits_%i.eps", nevent));
  canv->Print(Form("legendreHits_%i.png", nevent));
  delete canv;
 
  nevent++;
}

drawNode()

void drawNode

(

QuadTree *

node

)

const

Draw QuadTree node.

Definition at line 486 of file AxialHitQuadTreeProcessor.cc.

{
  std::vector<const CDCWireHit*> hits;
  for (Item* item : node->getItems()) {
    const CDCWireHit* wireHit = item->getPointer();
    hits.push_back(wireHit);
  }
  drawHits(hits);
}

fill() [1/2]

void fill ( const CandidateReceiver & candidateReceiver, int nHitsThreshold )

inlineinherited

Start filling the already created tree.

Parameters

candidateReceiver	the lambda function to call after a node was selected
nHitsThreshold	the threshold on the number of items

Definition at line 179 of file QuadTreeProcessor.h.

      {
        fill(candidateReceiver, nHitsThreshold, std::numeric_limits<AY>::max());
      }

fill() [2/2]

void fill ( const CandidateReceiver & candidateReceiver, int nHitsThreshold, float yLimit )

inlineinherited

Fill vector of QuadTree instances with hits.

Parameters

candidateReceiver	the lambda function to call after a node was selected
nHitsThreshold	the threshold on the number of items
yLimit	the threshold in the rho (curvature) variable

Definition at line 190 of file QuadTreeProcessor.h.

      {
        std::vector<QuadTree*> quadTrees = m_seededTrees;
        std::sort(quadTrees.begin(), quadTrees.end(), [](QuadTree * quadTree1, QuadTree * quadTree2) {
          return quadTree1->getNItems() > quadTree2->getNItems();
        });
 
        for (QuadTree* tree : quadTrees) {
          erase_remove_if(tree->getItems(), [](Item * hit) { return hit->isUsed(); });
          fillGivenTree(tree, candidateReceiver, nHitsThreshold, yLimit);
        }
      }

fillChildren()

void fillChildren ( QuadTree * node, const std::vector< Item * > & items )

inlineprivateinherited

This function is called by fillGivenTree and fills the items into the corresponding children.

For this the user-defined method isInNode is called.

Definition at line 281 of file QuadTreeProcessor.h.

      {
        // An item can be inserted into each child at most once
        const size_t neededSize = items.size();
        m_childPtrs.clear();
        m_childPtrs.reserve(node->getChildren().size());
        for (QuadTree& child : node->getChildren()) {
          child.reserveItems(neededSize);
          m_childPtrs.push_back(&child);
        }
 
        this->insertItemsInNodes(m_childPtrs, items);
 
        afterFillDebugHook(node->getChildren());
      }

fillGivenTree()

void fillGivenTree ( QuadTree * node, const CandidateReceiver & candidateReceiver, int nItemsThreshold, float yLimit )

inlineprivateinherited

Internal function to do the real quad tree search: fill the nodes, check which of the n*m bins we need to process further and go one level deeper.

Definition at line 208 of file QuadTreeProcessor.h.

      {
        if (node->getNItems() < nItemsThreshold) {
          return;
        }
 
        if ((node->getYMin() > yLimit) or (-node->getYMax() > yLimit)) {
          return;
        }
 
        if (isLeaf(node)) {
          callResultFunction(node, candidateReceiver);
          return;
        }
 
        if (node->getChildren().empty()) {
          this->createChildren(node, node->getChildren());
        }
 
        if (!node->checkFilled()) {
          fillChildren(node, node->getItems());
          node->setFilled();
        }
 
        // Kept on the stack - this function is called millions of times per event
        std::array<QuadTree*, c_maxNChildren> children;
        int nChildren = 0;
        for (QuadTree& child : node->getChildren()) {
          B2ASSERT("More children than the quad tree processor supports", nChildren < c_maxNChildren);
          children[nChildren++] = &child;
        }
        const auto compareNItems = [](const QuadTree * lhs, const QuadTree * rhs) {
          return lhs->getNItems() < rhs->getNItems();
        };
 
        // Explicitly count down the children
        for (int nRemaining = nChildren; nRemaining > 0; --nRemaining) {
          auto itHeaviestChild =
            std::max_element(children.begin(), children.begin() + nRemaining, compareNItems);
          QuadTree* heaviestChild = *itHeaviestChild;
          // Drop the heaviest child from the list keeping the order of the remaining ones
          std::move(itHeaviestChild + 1, children.begin() + nRemaining, itHeaviestChild);
          // After we have processed some children we need to get rid of the already used hits in all the children,
          // because this can change the number of items drastically
          erase_remove_if(heaviestChild->getItems(), [&](Item * hit) { return hit->isUsed(); });
          this->fillGivenTree(heaviestChild, candidateReceiver, nItemsThreshold, yLimit);
        }
      }

getAssignedItems()

std::vector< const TrackingUtilities::CDCWireHit * > getAssignedItems ( )

inlineinherited

Get items that have been assigned to the seed level The returned elements are unique even if items are assigned multiple times.

Definition at line 160 of file QuadTreeProcessor.h.

      {
        std::vector<const TrackingUtilities::CDCWireHit*> result;
        for (QuadTree* seededTree : m_seededTrees) {
          for (Item* item : seededTree->getItems()) {
            result.push_back(item->getPointer());
          }
        }
        std::sort(result.begin(), result.end());
        result.erase(std::unique(result.begin(), result.end()), result.end());
        return result;
      }

getCosSinLookupTable()

const LookupTable< Vector2D > & getCosSinLookupTable ( )

static

Get the standard lookup table containing equally spaces unit vectors (cos, sin)

It contains 2**16 + 1 sampling points between -pi and pi.

Definition at line 89 of file AxialHitQuadTreeProcessor.cc.

{
  static const int maxLevel = PrecisionUtil::getLookupGridLevel();
  static const int nBins = std::pow(2, maxLevel);
  static LookupTable<Vector2D> trigonometricLookUpTable(&Vector2D::Phi, nBins, -M_PI, M_PI);
  return trigonometricLookUpTable;
}

getDebugInformation()

const std::map< std::pair< long, float >, std::vector< Item * > > & getDebugInformation ( ) const

inlineinherited

Return the debug information if collected.

Definition at line 404 of file QuadTreeProcessor.h.

      {
        return m_debugOutputMap;
      }

getLastLevel()

int getLastLevel ( ) const

inlineprotectedinherited

Return the parameter last level.

Definition at line 382 of file QuadTreeProcessor.h.

      {
        return m_lastLevel;
      }

insertItemsInNodes() [1/2]

void insertItemsInNodes ( const std::vector< QuadTree * > & nodes, const std::vector< Item * > & items )

finalprotected

Insert the hits into the given nodes sharing the parts of the containment check that do not depend on the individual node.

The nodes handed to this function - the children of a common parent or the nodes of the seed level - always form a grid in the phase space, in the sense that their theta spans only take a few distinct values. All quantities that only depend on the hit and on the theta span of a node are therefore computed once per hit and distinct theta span instead of once for every hit and node pair.

Definition at line 262 of file AxialHitQuadTreeProcessor.cc.

{
  const size_t nNodes = nodes.size();
  if (nNodes == 0 or items.empty()) return;
 
  // Collect the geometry of the nodes and group them by their theta span.
  // The nodes handed here are either the four children of one node or the nodes of the seed
  // level, which both cover only a handful of distinct theta spans.
  m_thetaSpanCaches.clear();
  m_nodeCaches.resize(nNodes);
 
  for (size_t iNode = 0; iNode < nNodes; ++iNode) {
    QuadTree* node = nodes[iNode];
    NodeCache& nodeCache = m_nodeCaches[iNode];
    nodeCache.yMin = node->getYMin();
    nodeCache.yMax = node->getYMax();
    nodeCache.needsDerivativeCheck = node->getLevel() <= 4 and m_twoSidedPhaseSpace and
                                     nodeCache.yMin > -c_curlCurv and nodeCache.yMax < c_curlCurv;
 
    const long xMin = node->getXMin();
    const long xMax = node->getXMax();
 
    size_t iThetaSpan = 0;
    for (; iThetaSpan < m_thetaSpanCaches.size(); ++iThetaSpan) {
      if (m_thetaSpanCaches[iThetaSpan].xMin == xMin and m_thetaSpanCaches[iThetaSpan].xMax == xMax) break;
    }
    if (iThetaSpan == m_thetaSpanCaches.size()) {
      ThetaSpanCache thetaSpanCache;
      thetaSpanCache.xMin = xMin;
      thetaSpanCache.xMax = xMax;
      thetaSpanCache.thetaVecMin = &m_cosSinLookupTable->at(xMin);
      thetaSpanCache.thetaVecMax = &m_cosSinLookupTable->at(xMax);
      m_thetaSpanCaches.push_back(thetaSpanCache);
    }
    nodeCache.iThetaSpan = iThetaSpan;
  }
 
  const size_t nThetaSpans = m_thetaSpanCaches.size();
 
  for (Item* item : items) {
    if (item->isUsed()) continue;
 
    const CDCWireHit* wireHit = item->getPointer();
 
    // Quantities that only depend on the hit
    const double driftLength = wireHit->getRefDriftLength();
    const Vector2D pos2D = wireHit->getRefPos2D() - m_localOrigin;
    const double r2 = pos2D.normSquared() - driftLength * driftLength;
 
    // Quantities that only depend on the hit and the theta span of a node
    for (size_t iThetaSpan = 0; iThetaSpan < nThetaSpans; ++iThetaSpan) {
      ThetaSpanCache& thetaSpanCache = m_thetaSpanCaches[iThetaSpan];
      const Vector2D& thetaVecMin = *thetaSpanCache.thetaVecMin;
      const Vector2D& thetaVecMax = *thetaSpanCache.thetaVecMax;
 
      const float rHitMin = thetaVecMin.dot(pos2D);
      const float rHitMax = thetaVecMax.dot(pos2D);
 
      thetaSpanCache.rHitMinRight = rHitMin - driftLength;
      thetaSpanCache.rHitMaxRight = rHitMax - driftLength;
      thetaSpanCache.rHitMinLeft = rHitMin + driftLength;
      thetaSpanCache.rHitMaxLeft = rHitMax + driftLength;
 
      const float rHitMinExtr = thetaVecMin.cross(pos2D);
      const float rHitMaxExtr = thetaVecMax.cross(pos2D);
      thetaSpanCache.rHitMinExtr = rHitMinExtr;
      thetaSpanCache.rHitMaxExtr = rHitMaxExtr;
 
      // Same decision as checkDerivative()
      thetaSpanCache.derivativeOk = ((rHitMinExtr > 0) and (rHitMaxExtr * rHitMinExtr >= 0)) or
                                    (rHitMaxExtr * rHitMinExtr < 0);
 
      thetaSpanCache.hasExtremum = rHitMinExtr * rHitMaxExtr < 0.;
      if (thetaSpanCache.hasExtremum) {
        thetaSpanCache.extremumIsBetween = pos2D.isBetween(thetaVecMin, thetaVecMax);
      }
    }
 
    // Legendre curves at the extremum - only needed if some theta span contains the extremum
    bool extremumComputed = false;
    float rRight = 0;
    float rLeft = 0;
 
    for (size_t iNode = 0; iNode < nNodes; ++iNode) {
      const NodeCache& nodeCache = m_nodeCaches[iNode];
      const ThetaSpanCache& thetaSpanCache = m_thetaSpanCaches[nodeCache.iThetaSpan];
 
      // Check whether the hit lies in the forward direction
      if (nodeCache.needsDerivativeCheck and not thetaSpanCache.derivativeOk) continue;
 
      // get top and bottom borders of the node
      const float rMin = nodeCache.yMin * r2 / 2;
      const float rMax = nodeCache.yMax * r2 / 2;
 
      // Compare distance signs from the Legendre curves to the node
      // Check right
      if (not sameSign(rMin - thetaSpanCache.rHitMinRight,
                       rMin - thetaSpanCache.rHitMaxRight,
                       rMax - thetaSpanCache.rHitMinRight,
                       rMax - thetaSpanCache.rHitMaxRight)) {
        nodes[iNode]->insertItem(item);
        continue;
      }
 
      // Check left
      if (not sameSign(rMin - thetaSpanCache.rHitMinLeft,
                       rMin - thetaSpanCache.rHitMaxLeft,
                       rMax - thetaSpanCache.rHitMinLeft,
                       rMax - thetaSpanCache.rHitMaxLeft)) {
        nodes[iNode]->insertItem(item);
        continue;
      }
 
      // Check the extremum
      if (not thetaSpanCache.hasExtremum) continue;
      if (not thetaSpanCache.extremumIsBetween) continue;
 
      if (not extremumComputed) {
        const double r = pos2D.norm();
        rRight = r - driftLength;
        rLeft = r + driftLength;
        extremumComputed = true;
      }
 
      const bool crossesRight = (rMin - rRight) * (rMax - rRight) < 0;
      const bool crossesLeft = (rMin - rLeft) * (rMax - rLeft) < 0;
      if (crossesRight or crossesLeft) {
        nodes[iNode]->insertItem(item);
      }
    }
  }
}

insertItemsInNodes() [2/2]

virtual void insertItemsInNodes ( const std::vector< QuadTree * > & nodes, const std::vector< Item * > & items )

inlineprotectedvirtualinherited

Insert each of the given items into every one of the given nodes it belongs to.

The default implementation simply asks isInNode for each item and node pair. Processors may override it to share the parts of the containment check that only depend on the item or on a part of the node geometry over the given set of nodes. Any implementation has to insert the items into each node in the order in which they appear in the given item range.

Definition at line 351 of file QuadTreeProcessor.h.

      {
        for (Item* item : items) {
          if (item->isUsed()) continue;
 
          for (QuadTree* node : nodes) {
            if (isInNode(node, item->getPointer())) {
              node->insertItem(item);
            }
          }
        }
      }

isInNode() [1/2]

bool isInNode ( QuadTree * node, const TrackingUtilities::CDCWireHit * wireHit ) const

finalprotected

Check whether hit belongs to the quadtree node:

Parameters

node	quadtree node
wireHit	hit being checked

Returns

returns true if Legendre curve of the hit crosses (geometrically) borders of the node

Definition at line 194 of file AxialHitQuadTreeProcessor.cc.

{
  // Check whether the hit lies in the forward direction
  if (node->getLevel() <= 4 and m_twoSidedPhaseSpace and node->getYMin() > -c_curlCurv and
      node->getYMax() < c_curlCurv) {
    if (not checkDerivative(node, wireHit)) return false;
  }
 
  const double& driftLength = wireHit->getRefDriftLength();
  const Vector2D& pos2D = wireHit->getRefPos2D() - m_localOrigin;
  double r2 = pos2D.normSquared() - driftLength * driftLength;
 
  using Quadlet = std::array<std::array<float, 2>, 2>;
  Quadlet distRight{};
  Quadlet distLeft{};
 
  // get top and bottom borders of the node
  float rMin = node->getYMin() * r2 / 2;
  float rMax = node->getYMax() * r2 / 2;
 
  // get left and right borders of the node
  long thetaMin = node->getXMin();
  long thetaMax = node->getXMax();
 
  const Vector2D& thetaVecMin = m_cosSinLookupTable->at(thetaMin);
  const Vector2D& thetaVecMax = m_cosSinLookupTable->at(thetaMax);
 
  float rHitMin = thetaVecMin.dot(pos2D);
  float rHitMax = thetaVecMax.dot(pos2D);
 
  // compute Legendre curves at the left and right borders of the node
  float rHitMinRight = rHitMin - driftLength;
  float rHitMaxRight = rHitMax - driftLength;
 
  float rHitMinLeft = rHitMin + driftLength;
  float rHitMaxLeft = rHitMax + driftLength;
 
  // Compute distance from the Legendre curves to bottom and top borders of the node
  distRight[0][0] = rMin - rHitMinRight;
  distRight[0][1] = rMin - rHitMaxRight;
  distRight[1][0] = rMax - rHitMinRight;
  distRight[1][1] = rMax - rHitMaxRight;
 
  distLeft[0][0] = rMin - rHitMinLeft;
  distLeft[0][1] = rMin - rHitMaxLeft;
  distLeft[1][0] = rMax - rHitMinLeft;
  distLeft[1][1] = rMax - rHitMaxLeft;
 
  // Compare distance signs from Legendre curves to the node
  // Check right
  if (not sameSign(distRight[0][0], distRight[0][1], distRight[1][0], distRight[1][1])) {
    return true;
  }
 
  // Check left
  if (not sameSign(distLeft[0][0], distLeft[0][1], distLeft[1][0], distLeft[1][1])) {
    return true;
  }
 
  // Check the extremum
  float rHitMinExtr = thetaVecMin.cross(pos2D);
  float rHitMaxExtr = thetaVecMax.cross(pos2D);
  if (rHitMinExtr * rHitMaxExtr < 0.) return checkExtremum(node, wireHit);
 
  // Not contained
  return false;
}

isInNode() [2/2]

virtual bool isInNode ( QuadTree * node, const TrackingUtilities::CDCWireHit * item ) const

protectedpure virtualinherited

Implement that function if you want to provide a new processor.

It is called when filling the quad tree after creation. For every item in a node and every child node this function gets called and should decide, if the item should go into this child node or not.

Parameters

node	child node
item	item to be filled into the child node or not

Returns

true if this item belongs into this node.

isLeaf() [1/2]

bool isLeaf ( QuadTree * node ) const

finalprotected

lastLevel depends on curvature of the track candidate

Definition at line 121 of file AxialHitQuadTreeProcessor.cc.

{
  if (node->getLevel() <= 6) return false;
  if (node->getLevel() >= getLastLevel()) return true;
 
  const double nodeResolution = fabs(node->getYMin() - node->getYMax());
  const double meanCurv = (node->getYMax() + node->getYMin()) / 2;
 
  const double resolution = m_precisionFunction(meanCurv);
  if (resolution >= nodeResolution) return true;
 
  return false;
}

isLeaf() [2/2]

virtual bool isLeaf ( QuadTree * node ) const

inlineprotectedvirtualinherited

Function which checks if given node is leaf Implemented as virtual to keep possibility of changing lastLevel values depending on region is phase-space (i.e.

setting lastLevel as a function of Y-variable)

Definition at line 370 of file QuadTreeProcessor.h.

      {
        if (node->getLevel() >= m_lastLevel) {
          return true;
        } else {
          return false;
        }
      }

seed()

void seed ( const std::vector< const TrackingUtilities::CDCWireHit * > & datas )

inlineinherited

Fill in the items in the given vector.

They are transformed to QuadTreeItems internally.

Definition at line 113 of file QuadTreeProcessor.h.

      {
        // Create the items
        for (AData* data : datas) {
          m_items.emplace_back(data);
        }
 
        // Creating the seed level
        long nSeedBins = pow(2, m_seedLevel);
        m_seededTrees.reserve(nSeedBins * nSeedBins);
 
        // Expand the first levels to the seed sectors
        m_seededTrees.push_back(m_quadTree.get());
        std::vector<QuadTree*> nextSeededTrees;
 
        for (int level = 0; level < m_seedLevel; ++level) {
          for (QuadTree* node : m_seededTrees) {
            if (node->getChildren().empty()) {
              this->createChildren(node, node->getChildren());
            }
            for (QuadTree& child :  node->getChildren()) {
              nextSeededTrees.push_back(&child);
            }
          }
          std::swap(nextSeededTrees, m_seededTrees);
          nextSeededTrees.clear();
        }
 
        // Fill the seed level with the items
        m_itemPtrs.clear();
        m_itemPtrs.reserve(m_items.size());
        for (Item& item : m_items) {
          m_itemPtrs.push_back(&item);
        }
 
        for (QuadTree* seededTree : m_seededTrees) {
          seededTree->reserveItems(m_items.size());
        }
 
        this->insertItemsInNodes(m_seededTrees, m_itemPtrs);
      }

Member Data Documentation

c_curlCurv

const double c_curlCurv = 0.02

private

The curvature above which the trajectory is considered a curler.

Definition at line 202 of file AxialHitQuadTreeProcessor.h.

c_maxNChildren

const int c_maxNChildren

staticprivateinherited

Upper bound on the number of children of a node - they are kept on the stack.

Definition at line 47 of file QuadTreeProcessor.h.

m_childPtrs

std::vector<QuadTree*> m_childPtrs

protectedinherited

Reusable buffer with pointers to the children of the node currently being filled.

Definition at line 420 of file QuadTreeProcessor.h.

m_cosSinLookupTable

const TrackingUtilities::LookupTable<TrackingUtilities::Vector2D>* m_cosSinLookupTable

private

Pinned lookup table for precomputed cosine and sine values.

Definition at line 199 of file AxialHitQuadTreeProcessor.h.

m_debugOutput

bool m_debugOutput

privateinherited

A flag to control the creation of the debug output.

Definition at line 437 of file QuadTreeProcessor.h.

m_debugOutputMap

std::map<std::pair<long, float>, std::vector<Item*> > m_debugOutputMap

privateinherited

The calculated debug map.

Definition at line 440 of file QuadTreeProcessor.h.

m_itemPtrs

std::vector<Item*> m_itemPtrs

protectedinherited

Reusable buffer with pointers to all items - used to seed the tree.

Definition at line 417 of file QuadTreeProcessor.h.

m_items

std::deque<Item> m_items

protectedinherited

Storage space for the items that are referenced by the quad tree nodes.

Definition at line 414 of file QuadTreeProcessor.h.

m_lastLevel

int m_lastLevel

privateinherited

The last level to be filled.

Definition at line 431 of file QuadTreeProcessor.h.

m_localOrigin

TrackingUtilities::Vector2D m_localOrigin

private

Local origin on which the phase space coordinates are centered.

Definition at line 196 of file AxialHitQuadTreeProcessor.h.

m_nodeCaches

std::vector<NodeCache> m_nodeCaches

private

Reusable buffer with the per node quantities - one entry per node.

Definition at line 189 of file AxialHitQuadTreeProcessor.h.

m_precisionFunction

PrecisionUtil::PrecisionFunction m_precisionFunction

private

Lambda which holds resolution function for the quadtree.

Definition at line 193 of file AxialHitQuadTreeProcessor.h.

m_quadTree

std::unique_ptr<QuadTree> m_quadTree

protectedinherited

The quad tree we work with.

Definition at line 411 of file QuadTreeProcessor.h.

m_seededTrees

std::vector<QuadTree*> m_seededTrees

protectedinherited

Vector of QuadTrees QuadTree instances (which are filled in the vector) cover the whole Legendre phase-space; each instance is processes independently.

Definition at line 427 of file QuadTreeProcessor.h.

m_seedLevel

int m_seedLevel

privateinherited

The first level to be filled, effectively skip forward to this higher granularity level.

Definition at line 434 of file QuadTreeProcessor.h.

m_thetaSpanCaches

std::vector<ThetaSpanCache> m_thetaSpanCaches

private

Reusable buffer with the per theta span quantities - one entry per distinct theta span.

Definition at line 186 of file AxialHitQuadTreeProcessor.h.

m_twoSidedPhaseSpace

bool m_twoSidedPhaseSpace

private

Indicator whether the two sided phases space insertion check should be used This option should automatically split back to back tracks in the low curvature regions.

Definition at line 208 of file AxialHitQuadTreeProcessor.h.

---

The documentation for this class was generated from the following files:

• tracking/trackFindingCDC/legendre/quadtree/include/AxialHitQuadTreeProcessor.h
• tracking/trackFindingCDC/legendre/quadtree/src/AxialHitQuadTreeProcessor.cc