Class to represent the CDC NDFinder. More...

#include <NDFinder.h>

Classes
struct	ContributionInfo
	Collection of the hit contribution information needed for the hit to cluster relations. More...

struct	SectorBinning
	Data type to collect a binning. More...

struct	WireInfo
	Collection of the hit information needed for the hit representations. More...

Public Member Functions
void	init (const NDFinderParameters &ndFinderParameters)
	Initialization of the NDFinder (parameters and lookup tables)

void	reset ()
	Reset the NDFinder data structure to process next event.

void	addHit (const HitInfo &hitInfo)
	Add the hit info of a single track segment to the NDFinder.

void	findTracks ()
	Main function for track finding.

std::vector< NDFinderTrack > *	getFinderTracks ()
	Retreive the results.

Protected Member Functions
void	initLookUpArrays ()
	Initialize the arrays LUT arrays.

void	initHitToSectorMap ()
	Fills the m_hitToSectorIDs (see below) array with the hit to orientation/sectorWire/sectorID relations.

void	loadCompressedHitReps (const std::string &fileName, const SectorBinning &compBins, c5array &compHitsToWeights)
	Fills the m_compAxialHitReps/m_compStereoHitReps (see below) arrays with the hit representations (hits to weights)

void	fillExpandedHitReps (const SectorBinning &compBins, const c5array &compHitsToWeights, c5array &expHitsToWeights)
	Fills the m_expAxialHitReps/m_expStereoHitReps (see below) arrays with the expanded hit representations (hits to weights)

void	processHitForHoughSpace (const unsigned short hitIdx)
	Process a single axial or stereo hit for the Hough space.

void	writeHitToHoughSpace (const WireInfo &hitInfo, const c5array &expHitsToWeights)
	Write (add) a single hit (Hough curve) to the Hough space.

void	runTrackFinding ()
	Core track finding logic in the constructed Hough space.

std::vector< SimpleCluster >	relateHitsToClusters (std::vector< SimpleCluster > &clusters)
	Relate the hits in the peak of the cluster to the cluster. Applies a cut on the clusters.

std::vector< std::vector< unsigned short > >	getHitsVsClustersTable (const std::vector< SimpleCluster > &clusters)
	Create hits to clusters confusion matrix.

unsigned short	getHitContribution (const cell_index &peakCell, const unsigned short hitIdx)
	Returns the hit contribution of a TS at a certain cluster cell (= peak/maximum cell)

std::vector< ContributionInfo >	extractContributionInfos (const std::vector< unsigned short > &clusterHits)
	Extract relevant hit information (hitIdx, contribution, super layer, drift time)

int	getMaximumHitInSuperLayer (const std::vector< ContributionInfo > &contributionInfos, unsigned short superLayer)
	Find the hit with the maximum contribution in a given super layer.

bool	checkHitSuperLayers (const SimpleCluster &cluster)
	Cut on the number of hit axial/stereo super layers.

std::array< double, 3 >	calculateCenterOfGravity (const SimpleCluster &cluster)
	Calculate the center of gravity (weighted mean) for the track parameters.

std::array< double, 3 >	getTrackParameterEstimate (const std::array< double, 3 > &centerOfGravity)
	Transform the center of gravity (cells) into the estimated track parameters.

std::array< double, 3 >	transformTrackParameters (const std::array< double, 3 > &estimatedParameters)
	Transform to physical units.

Static Protected Member Functions
static double	getTrackRadius (double transverseMomentum)
	Transverse momentum (which is 1/omega, in GeV/c) to radius (in cm)

Private Attributes
std::vector< NDFinderTrack >	m_ndFinderTracks
	Result: Vector of the found tracks.

std::vector< unsigned short >	m_hitIDs
	TS-IDs of the hits in the current event: Elements = [0,2335] for 2336 TS in total.

std::vector< unsigned short >	m_hitSLIDs
	SL-IDs of the hits in the current event: Elements = Super layer number in [0,1,...,8].

std::vector< unsigned short >	m_priorityWirePos
	Priority positon within the TS. Elements basf2: [0,3] first, left, right, no hit.

std::vector< short >	m_priorityWireTime
	Drift time of the priority wire.

unsigned short	m_nHits {0}
	Counter for the number of hits in the current event.

NDFinderParameters	m_ndFinderParams
	Configuration parameters of the 3DFinder.

Clusterizend	m_clusterer
	Clustering module.

c2array *	m_hitToSectorIDs = nullptr
	m_hitToSectorIDs: 2D array mapping TS-ID ([0, 2335]) to:

c5array *	m_compAxialHitReps = nullptr
	m_compAxialHitReps/m_compStereoHitReps (~ Compressed in phi (width, start, values)) 5D array mapping:

c5array *	m_compStereoHitReps = nullptr

c5array *	m_expAxialHitReps = nullptr
	m_expAxialHitReps/m_expStereoHitReps (~ expansion of the compressed representations) 5D array mapping:

c5array *	m_expStereoHitReps = nullptr

c3array *	m_houghSpace = nullptr
	The complete Hough space with the size [m_nOmega, m_nPhi, m_nCot].

Static Private Attributes
static constexpr unsigned short	m_nTS = 2336
	Number of track segments.

static constexpr unsigned short	m_nSL = 9
	Number of super layers.

static constexpr unsigned short	m_nAxial = 41
	Number of unique axial track segments.

static constexpr unsigned short	m_nStereo = 32
	Number of unique stereo track segments.

static constexpr unsigned short	m_nPrio = 3
	Number of priority wires.

static constexpr unsigned short	m_nOmega = 40
	Bins in the phi dimension.

static constexpr unsigned short	m_nPhi = 384
	Bins in the omega dimension.

static constexpr unsigned short	m_nCot = 9
	Bins in the cot dimension.

static constexpr unsigned short	m_phiGeo = 32
	Repetition of the wire pattern.

static constexpr unsigned short	m_nPhiSector = m_nPhi / m_phiGeo
	Bins of one phi sector (12)

static constexpr unsigned short	m_nPhiComp = 15
	Bins of compressed phi: phi_start, phi_width, phi_0, ..., phi_12.

static constexpr SectorBinning	m_compAxialBins = {m_nOmega, m_nPhiComp, 1, m_nAxial, m_nPrio}
	40, 15, 1, 41, 3

static constexpr SectorBinning	m_compStereoBins = {m_nOmega, m_nPhiComp, m_nCot, m_nStereo, m_nPrio}
	40, 15, 9, 32, 3

static constexpr std::array< double, 2 >	m_omegaRange = {-4., 4.}
	1/4 = 0.25 GeV (minimum transverse momentum)

static constexpr std::array< double, 2 >	m_phiRange = {0., 11.25}
	One phi sector (360/32)

static constexpr std::array< double, 2 >	m_cotRange = {2.3849627654510415, -1.0061730449796316}
	=> theta in [22.75, 135.18]

static constexpr double	m_binSizeOmega = (m_omegaRange[1] - m_omegaRange[0]) / m_nOmega
	0.2

static constexpr double	m_binSizePhi = (m_phiRange[1] - m_phiRange[0]) / m_nPhiSector
	0.9375

static constexpr double	m_binSizeCot = (m_cotRange[1] - m_cotRange[0]) / m_nCot
	-0.377

static constexpr std::array< std::array< double, 2 >, 3 >	m_acceptanceRanges = {m_omegaRange, m_phiRange, m_cotRange}

static constexpr std::array< double, 3 >	m_binSizes = {m_binSizeOmega, m_binSizePhi, m_binSizeCot}

Detailed Description

Class to represent the CDC NDFinder.

Definition at line 99 of file NDFinder.h.

Constructor & Destructor Documentation

◆ ~NDFinder()

virtual ~NDFinder ( )

inlinevirtual

Definition at line 129 of file NDFinder.h.

    {
      delete m_hitToSectorIDs;
      delete m_compAxialHitReps;
      delete m_compStereoHitReps;
      delete m_expAxialHitReps;
      delete m_expStereoHitReps;
      delete m_houghSpace;
    }

Member Function Documentation

◆ addHit()

void addHit ( const HitInfo & hitInfo )

Add the hit info of a single track segment to the NDFinder.

Definition at line 59 of file NDFinder.cc.

{
  // Priority position from the track segment finder: 0 = no hit, 3 = 1st priority, 1 = 2nd right, 2 = 2nd left
  if (hitInfo.hitPrioPos > 0) { // Skip "no hit" case
    m_hitIDs.push_back(hitInfo.hitID);
    m_hitSLIDs.push_back(hitInfo.hitSLID);
    m_priorityWirePos.push_back(3 - hitInfo.hitPrioPos); // 0 = 1st priority, 1 = 2nd left, 2 = 2nd right
    m_priorityWireTime.push_back(hitInfo.hitPrioTime);
    ++m_nHits;
  }
}

◆ calculateCenterOfGravity()

std::array< double, 3 > calculateCenterOfGravity ( const SimpleCluster & cluster )

protected

Calculate the center of gravity (weighted mean) for the track parameters.

Definition at line 430 of file NDFinder.cc.

{
  double weightedSumOmega = 0.;
  double weightedSumPhi = 0.;
  double weightedSumCot = 0.;
  unsigned int weightSum = 0;
  std::vector<cell_index> clusterCells = cluster.getCells();
  const int peakPhi = clusterCells.front()[1];
  for (const cell_index& cellIdx : clusterCells) {
    const unsigned short cellWeight = (*m_houghSpace)[cellIdx[0]][cellIdx[1]][cellIdx[2]];
    weightedSumOmega += cellIdx[0] * cellWeight;
    weightedSumCot += cellIdx[2] * cellWeight;
    // Calculate relative phi position to peak
    int phi = cellIdx[1];
    int delta = phi - peakPhi;
    // Wrap-around correction to stay within [-m_nPhi/2, +m_nPhi/2]
    if (delta > m_nPhi / 2) {
      delta -= m_nPhi;
    } else if (delta < -m_nPhi / 2) {
      delta += m_nPhi;
    }
    double unwrappedPhi = peakPhi + delta;
    weightedSumPhi += unwrappedPhi * cellWeight;
    weightSum += cellWeight;
  }
  weightedSumOmega /= weightSum;
  weightedSumCot /= weightSum;
  weightedSumPhi /= weightSum;
  if (weightedSumPhi >= m_nPhi) weightedSumPhi -= m_nPhi;
  if (weightedSumPhi < 0) weightedSumPhi += m_nPhi;
  std::array<double, 3> centerOfGravity = {weightedSumOmega, weightedSumPhi, weightedSumCot};
  return centerOfGravity;
}

◆ checkHitSuperLayers()

bool checkHitSuperLayers ( const SimpleCluster & cluster )

protected

Cut on the number of hit axial/stereo super layers.

Definition at line 410 of file NDFinder.cc.

{
  std::vector<unsigned short> clusterHits = cluster.getClusterHits();
  std::set<unsigned short> uniqueSLNumbers;
  // Add all hit super layers
  for (unsigned short hitIdx : clusterHits) {
    uniqueSLNumbers.insert(m_hitSLIDs[hitIdx]);
  }
  // Calculate the number of hit axial and stereo super layers
  unsigned short nSL = uniqueSLNumbers.size();
  uniqueSLNumbers.insert({0, 2, 4, 6, 8});
  unsigned short withAxialSLs = uniqueSLNumbers.size();
  unsigned short axialNumber = 5 - (withAxialSLs - nSL);
  unsigned short stereoNumber = nSL - axialNumber;
  // Cut away all clusters that do have enough hit super layers
  bool isValid = axialNumber >= m_ndFinderParams.minSuperAxial && stereoNumber >= m_ndFinderParams.minSuperStereo;
  return isValid;
}

◆ extractContributionInfos()

std::vector< NDFinder::ContributionInfo > extractContributionInfos ( const std::vector< unsigned short > & clusterHits )

protected

Extract relevant hit information (hitIdx, contribution, super layer, drift time)

Definition at line 363 of file NDFinder.cc.

{
  std::vector<ContributionInfo> contributionInfos;
  for (unsigned short hitIdx = 0; hitIdx < m_hitIDs.size(); ++hitIdx) {
    unsigned short contribution = clusterHits[hitIdx];
    if (contribution > 0) {
      ContributionInfo contributionInfo = {
        hitIdx,
        contribution,
        m_hitSLIDs[hitIdx],
        m_priorityWireTime[hitIdx]
      };
      contributionInfos.push_back(contributionInfo);
    }
  }
  return contributionInfos;
}

◆ fillExpandedHitReps()

void fillExpandedHitReps	(	const SectorBinning &	compBins,
		const c5array &	compHitsToWeights,
		c5array &	expHitsToWeights )

protected

Fills the m_expAxialHitReps/m_expStereoHitReps (see below) arrays with the expanded hit representations (hits to weights)

Definition at line 159 of file NDFinder.cc.

{
  for (c5index hitID = 0; hitID < compBins.nHitIDs; ++hitID) {
    for (c5index priorityWire = 0; priorityWire < compBins.nPriorityWires; ++priorityWire) {
      for (c5index omegaIdx = 0; omegaIdx < compBins.omega; ++omegaIdx) {
        for (c5index cotIdx = 0; cotIdx < compBins.cot; ++cotIdx) {
          unsigned short phiStart = compHitsToWeights[hitID][priorityWire][omegaIdx][0][cotIdx];
          unsigned short nPhiEntries = compHitsToWeights[hitID][priorityWire][omegaIdx][1][cotIdx];
          for (c5index phiEntry = 0; phiEntry < nPhiEntries; ++phiEntry) {
            unsigned short houghPhiIdx = (phiStart + phiEntry) % m_nPhi; // houghPhiIdx goes now over the complete Hough space
            expHitsToWeights[hitID][priorityWire][omegaIdx][houghPhiIdx][cotIdx] =
              compHitsToWeights[hitID][priorityWire][omegaIdx][phiEntry + 2][cotIdx];
            if (compBins.cot == 1) { // Axial case: expand the same curve in all cot bins
              for (c5index axialCotIdx = 1; axialCotIdx < m_nCot; ++axialCotIdx) {
                expHitsToWeights[hitID][priorityWire][omegaIdx][houghPhiIdx][axialCotIdx] =
                  compHitsToWeights[hitID][priorityWire][omegaIdx][phiEntry + 2][cotIdx];
              }
            }
          }
        }
      }
    }
  }
}

◆ findTracks()

void findTracks ( )

Main function for track finding.

Definition at line 204 of file NDFinder.cc.

{
  // Build the Hough plane by summing up all single hit contributions
  for (unsigned short hitIdx = 0; hitIdx < m_nHits; ++hitIdx) {
    processHitForHoughSpace(hitIdx);
  }
  runTrackFinding();
}

◆ getFinderTracks()

std::vector< NDFinderTrack > * getFinderTracks ( )

inline

Retreive the results.

Definition at line 148 of file NDFinder.h.

148{ return &m_ndFinderTracks; }

◆ getHitContribution()

unsigned short getHitContribution	(	const cell_index &	peakCell,
		const unsigned short	hitIdx )

protected

Returns the hit contribution of a TS at a certain cluster cell (= peak/maximum cell)

Definition at line 336 of file NDFinder.cc.

{
  unsigned short contribution = 0;
  unsigned short omegaIdx = peakCell[0];
  unsigned short houghPhiIdx = peakCell[1];
  unsigned short cotIdx = peakCell[2];
 
  const c2array& hitToSectorIDs = *m_hitToSectorIDs;
  unsigned short orientation = hitToSectorIDs[m_hitIDs[hitIdx]][0];
  unsigned short relativeWireID = hitToSectorIDs[m_hitIDs[hitIdx]][1];
  unsigned short relativeSectorID = hitToSectorIDs[m_hitIDs[hitIdx]][2];
  unsigned short phiSectorStart = relativeSectorID * m_nPhiSector;
  unsigned short priorityWire = m_priorityWirePos[hitIdx];
 
  // Inverse Hough transformation (inverse of writeHitToHoughSpace method)
  unsigned short phiIdx = (houghPhiIdx - phiSectorStart + m_nPhi) % m_nPhi;
 
  if (orientation == 1) { // Axial TS
    contribution = (*m_expAxialHitReps)[relativeWireID][priorityWire][omegaIdx][phiIdx][cotIdx];
  } else { // Stereo TS
    contribution = (*m_expStereoHitReps)[relativeWireID][priorityWire][omegaIdx][phiIdx][cotIdx];
  }
 
  return contribution;
}

◆ getHitsVsClustersTable()

std::vector< std::vector< unsigned short > > getHitsVsClustersTable ( const std::vector< SimpleCluster > & clusters )

protected

Create hits to clusters confusion matrix.

Definition at line 318 of file NDFinder.cc.

{
  // Creates a (Number Clusters)x(Number Hits) matrix
  std::vector<unsigned short> hitElem(m_nHits, 0);
  std::vector<std::vector<unsigned short>> hitsVsClusters(clusters.size(), hitElem);
  // Fill the matrix with all the hit contributions
  for (unsigned short clusterIdx = 0; clusterIdx < clusters.size(); ++clusterIdx) {
    SimpleCluster cluster = clusters[clusterIdx];
    cell_index peakCell = cluster.getPeakCell();
    for (unsigned short hitIdx = 0; hitIdx < m_hitIDs.size(); ++hitIdx) {
      unsigned short contribution = getHitContribution(peakCell, hitIdx);
      hitsVsClusters[clusterIdx][hitIdx] = contribution;
    }
  }
  return hitsVsClusters;
}

◆ getMaximumHitInSuperLayer()

int getMaximumHitInSuperLayer	(	const std::vector< ContributionInfo > &	contributionInfos,
		unsigned short	superLayer )

protected

Find the hit with the maximum contribution in a given super layer.

Definition at line 382 of file NDFinder.cc.

{
  std::vector<std::vector<int>> contributionsInSL;
  for (const ContributionInfo& contributionInfo : contributionInfos) {
    if (contributionInfo.superLayer == superLayer) {
      unsigned short hitIdx = contributionInfo.hitIndex;
      unsigned short contribution = contributionInfo.contribution;
      short priorityTime = contributionInfo.priorityTime;
      contributionsInSL.push_back({hitIdx, contribution, priorityTime});
    }
  }
  if (contributionsInSL.empty()) return -1;
  // Sort by drift time
  std::sort(contributionsInSL.begin(), contributionsInSL.end(),
  [](const std::vector<int>& a, const std::vector<int>& b) { return a[2] < b[2]; });
  // Find max contribution
  int maximumHit = contributionsInSL[0][0];
  int maximumContribution = contributionsInSL[0][1];
  for (const auto& hit : contributionsInSL) {
    if (hit[1] > maximumContribution) {
      maximumHit = hit[0];
      maximumContribution = hit[1];
    }
  }
  return maximumHit;
}

◆ getTrackParameterEstimate()

std::array< double, 3 > getTrackParameterEstimate ( const std::array< double, 3 > & centerOfGravity )

protected

Transform the center of gravity (cells) into the estimated track parameters.

Definition at line 465 of file NDFinder.cc.

{
  std::array<double, 3> estimatedParameters;
  for (unsigned short dimension = 0; dimension < 3; ++dimension) {
    double trackParameter = m_acceptanceRanges[dimension][0] + (centerOfGravity[dimension] + 0.5) * m_binSizes[dimension];
    estimatedParameters[dimension] = trackParameter;
  }
  return transformTrackParameters(estimatedParameters);
}

◆ getTrackRadius()

static double getTrackRadius ( double transverseMomentum )

inlinestaticprotected

Transverse momentum (which is 1/omega, in GeV/c) to radius (in cm)

Definition at line 185 of file NDFinder.h.

185{ return transverseMomentum * 1e11 / (3e8 * 1.5); }

◆ init()

void init ( const NDFinderParameters & ndFinderParameters )

Initialization of the NDFinder (parameters and lookup tables)

Definition at line 27 of file NDFinder.cc.

{
  m_ndFinderParams = ndFinderParameters;
 
  // Initialization of the pointer arrays, fills hit to sector LUT
  initLookUpArrays();
  initHitToSectorMap();
 
  // Load the axial and stereo track to hit relations from file.
  loadCompressedHitReps(m_ndFinderParams.axialFile, m_compAxialBins, *m_compAxialHitReps);
  loadCompressedHitReps(m_ndFinderParams.stereoFile, m_compStereoBins, *m_compStereoHitReps);
 
  // Fills the expanded hit representations (from compressed hits to weights)
  fillExpandedHitReps(m_compAxialBins, *m_compAxialHitReps, *m_expAxialHitReps);
  fillExpandedHitReps(m_compStereoBins, *m_compStereoHitReps, *m_expStereoHitReps);
 
  // Reset the NDFinder data structure to process next event
  reset();
 
  // Parameters necessary for the clustering algorithm
  ClustererParameters clustererParams = {
    ndFinderParameters.iterations,
    ndFinderParameters.omegaTrim,
    ndFinderParameters.phiTrim,
    m_nOmega,
    m_nPhi,
    m_nCot
  };
  m_clusterer = Clusterizend(clustererParams);
}

◆ initHitToSectorMap()

void initHitToSectorMap ( )

protected

Fills the m_hitToSectorIDs (see below) array with the hit to orientation/sectorWire/sectorID relations.

Definition at line 91 of file NDFinder.cc.

{
  // Number of first priority wires in each super layer (TS per SL)
  constexpr std::array<unsigned short, 9> nWires = {160, 160, 192, 224, 256, 288, 320, 352, 384};
  // Number of priority wires (= number of TS) per SL in a single (1/32) phi sector
  std::vector<unsigned short> wiresPerSector;
  // Lookup table: Maps the TS id to the SL number
  std::vector<unsigned short> hitToSuperLayer;
  // Integrated number of priority wires for each SL
  std::vector<unsigned short> cumulativeWires = {0};
  // Integrated number of sector priority wires for each SL (Axial even, Stereo odd)
  std::vector<unsigned short> cumulativeSectorWires = {0, 0};
  for (unsigned short superLayer = 0; superLayer < m_nSL; ++superLayer) {
    wiresPerSector.push_back(nWires[superLayer] / m_phiGeo);
    for (unsigned short _ = 0; _ < nWires[superLayer]; ++_) {
      hitToSuperLayer.push_back(superLayer);
    }
    cumulativeWires.push_back(cumulativeWires[superLayer] + nWires[superLayer]);
    cumulativeSectorWires.push_back(cumulativeSectorWires[superLayer] + nWires[superLayer] / m_phiGeo);
  }
  for (unsigned short hit = 0; hit < m_nTS; ++hit) {
    unsigned short superLayer = hitToSuperLayer[hit];
    bool isAxial = (superLayer % 2 == 0);
    unsigned short wireIDinSL = hit - cumulativeWires[superLayer];
    unsigned short wireIDinSector = wireIDinSL % wiresPerSector[superLayer];
    unsigned short relativeWireIDinSector = cumulativeSectorWires[superLayer] + wireIDinSector;
    unsigned short relativeSectorIDinSuperLayer = static_cast<unsigned short>(floor(wireIDinSL / wiresPerSector[superLayer]));
    c2array& hitToSectorIDs = *m_hitToSectorIDs;
    hitToSectorIDs[hit][0] = static_cast<unsigned short>(isAxial);
    hitToSectorIDs[hit][1] = relativeWireIDinSector;
    hitToSectorIDs[hit][2] = relativeSectorIDinSuperLayer;
  }
}

◆ initLookUpArrays()

void initLookUpArrays ( )

protected

Initialize the arrays LUT arrays.

Definition at line 72 of file NDFinder.cc.

{
  // Shapes of the arrays holding the hit patterns
  const std::array<c2index, 2> shapeHitToSectorIDs = {{m_nTS, m_nPrio}};
  const std::array<c5index, 5> shapeCompAxialHitReps = {{m_nAxial, m_nPrio, m_nOmega, m_nPhiComp, 1}};
  const std::array<c5index, 5> shapeCompStereoHitReps = {{m_nStereo, m_nPrio, m_nOmega, m_nPhiComp, m_nCot}};
  const std::array<c5index, 5> shapeExpAxialHitReps = {{m_nAxial, m_nPrio, m_nOmega, m_nPhi, m_nCot}};
  const std::array<c5index, 5> shapeExpStereoHitReps = {{m_nStereo, m_nPrio, m_nOmega, m_nPhi, m_nCot}};
  const std::array<c3index, 3> shapeHough = {{m_nOmega, m_nPhi, m_nCot}};
 
  m_hitToSectorIDs    = new c2array(shapeHitToSectorIDs);
  m_compAxialHitReps  = new c5array(shapeCompAxialHitReps);
  m_compStereoHitReps = new c5array(shapeCompStereoHitReps);
  m_expAxialHitReps   = new c5array(shapeExpAxialHitReps);
  m_expStereoHitReps  = new c5array(shapeExpStereoHitReps);
  m_houghSpace        = new c3array(shapeHough);
}

◆ loadCompressedHitReps()

void loadCompressedHitReps	(	const std::string &	fileName,
		const SectorBinning &	compBins,
		c5array &	compHitsToWeights )

protected

Fills the m_compAxialHitReps/m_compStereoHitReps (see below) arrays with the hit representations (hits to weights)

Definition at line 126 of file NDFinder.cc.

{
  // Array of the entries in trg/cdc/data/ndFinderArray*Comp.txt.gz
  std::vector<unsigned short> flatArray;
  std::ifstream arrayFileGZ(fileName, std::ios_base::in | std::ios_base::binary);
  if (!arrayFileGZ.is_open()) {
    B2ERROR("Could not open array file: " << fileName);
    return;
  }
  boost::iostreams::filtering_istream arrayStream;
  arrayStream.push(boost::iostreams::gzip_decompressor());
  arrayStream.push(arrayFileGZ);
  unsigned short uline;
  while (arrayStream >> uline) {
    flatArray.push_back(uline);
  }
  arrayFileGZ.close();
  unsigned long arrayIndex = 0;
  for (c5index hitID = 0; hitID < compBins.nHitIDs; ++hitID) {
    for (c5index priorityWire = 0; priorityWire < compBins.nPriorityWires; ++priorityWire) {
      for (c5index omegaIdx = 0; omegaIdx < compBins.omega; ++omegaIdx) {
        for (c5index phiIdx = 0; phiIdx < compBins.phi; ++phiIdx) {
          for (c5index cotIdx = 0; cotIdx < compBins.cot; ++cotIdx) {
            compHitsToWeights[hitID][priorityWire][omegaIdx][phiIdx][cotIdx] = flatArray[arrayIndex];
            ++arrayIndex;
          }
        }
      }
    }
  }
}

◆ processHitForHoughSpace()

void processHitForHoughSpace ( const unsigned short hitIdx )

protected

Process a single axial or stereo hit for the Hough space.

Definition at line 214 of file NDFinder.cc.

{
  const c2array& hitToSectorIDs = *m_hitToSectorIDs;
  unsigned short orientation = hitToSectorIDs[m_hitIDs[hitIdx]][0];
  unsigned short relativeWireID = hitToSectorIDs[m_hitIDs[hitIdx]][1];
  unsigned short relativeSectorID = hitToSectorIDs[m_hitIDs[hitIdx]][2];
  unsigned short phiSectorStart = relativeSectorID * m_nPhiSector;
  unsigned short priorityWire = m_priorityWirePos[hitIdx];
 
  WireInfo wireInfo = {relativeWireID, phiSectorStart, priorityWire};
  if (orientation == 1) {
    writeHitToHoughSpace(wireInfo, *m_expAxialHitReps);
  } else {
    writeHitToHoughSpace(wireInfo, *m_expStereoHitReps);
  }
}

◆ relateHitsToClusters()

std::vector< SimpleCluster > relateHitsToClusters ( std::vector< SimpleCluster > & clusters )

protected

Relate the hits in the peak of the cluster to the cluster. Applies a cut on the clusters.

Definition at line 296 of file NDFinder.cc.

{
  std::vector<std::vector<unsigned short>> hitsVsClusters = getHitsVsClustersTable(clusters);
  std::vector<SimpleCluster> goodClusters;
  if (hitsVsClusters.empty()) return goodClusters;
 
  for (unsigned short clusterIdx = 0; clusterIdx < hitsVsClusters.size(); ++clusterIdx) {
    std::vector<ContributionInfo> contributionInfos = extractContributionInfos(hitsVsClusters[clusterIdx]);
    for (unsigned short superLayer = 0; superLayer < 9; ++superLayer) {
      int maximumHit = getMaximumHitInSuperLayer(contributionInfos, superLayer);
      if (maximumHit >= 0) { // there exists a hit
        clusters[clusterIdx].addHitToCluster(static_cast<unsigned short>(maximumHit));
      }
    }
    if (checkHitSuperLayers(clusters[clusterIdx])) {
      goodClusters.push_back(clusters[clusterIdx]);
    }
  }
  return goodClusters;
}

◆ reset()

void reset ( )

Reset the NDFinder data structure to process next event.

Definition at line 185 of file NDFinder.cc.

{
  // Clear the vector of the found tracks
  m_ndFinderTracks.clear();
 
  // Clear the hit informations
  m_hitIDs.clear();
  m_hitSLIDs.clear();
  m_nHits = 0;
  m_priorityWirePos.clear();
  m_priorityWireTime.clear();
 
  // Create a new Hough space
  delete m_houghSpace;
  std::array<c3index, 3> shapeHough = {{m_nOmega, m_nPhi, m_nCot}};
  m_houghSpace = new c3array(shapeHough);
}

◆ runTrackFinding()

void runTrackFinding ( )

protected

Core track finding logic in the constructed Hough space.

Definition at line 247 of file NDFinder.cc.

{
  c3array& houghSpace = *m_houghSpace;
  m_clusterer.setNewPlane(houghSpace);
  std::vector<SimpleCluster> allClusters = m_clusterer.makeClusters();
  std::vector<SimpleCluster> validClusters = relateHitsToClusters(allClusters);
 
  for (SimpleCluster& cluster : validClusters) {
    std::array<double, 3> centerOfGravity = calculateCenterOfGravity(cluster);
    std::array<double, 3> estimatedParameters = getTrackParameterEstimate(centerOfGravity);
 
    // Readouts for the 3DFinderInfo class for analysis (Hough space + cluster info)
    std::vector<ROOT::Math::XYZVector> readoutHoughSpace;
    std::vector<ROOT::Math::XYZVector> readoutCluster;
 
    if (m_ndFinderParams.storeAdditionalReadout) {
      // Readout of the complete Hough space
      for (c3index omegaIdx = 0; omegaIdx < m_nOmega; ++omegaIdx) {
        for (c3index phiIdx = 0; phiIdx < m_nPhi; ++phiIdx) {
          for (c3index cotIdx = 0; cotIdx < m_nCot; ++cotIdx) {
            unsigned short element = houghSpace[omegaIdx][phiIdx][cotIdx];
            readoutHoughSpace.push_back(ROOT::Math::XYZVector(element, 0, 0));
          }
        }
      }
      // Readout of the peak cluster weight
      unsigned int peakWeight = cluster.getPeakWeight();
      readoutCluster.push_back(ROOT::Math::XYZVector(peakWeight, 0, 0));
      // Readout of the number of cluster cells
      unsigned short nCells = cluster.getCells().size();
      readoutCluster.push_back(ROOT::Math::XYZVector(nCells, 0, 0));
      // Readout of the cluster center of gravity
      readoutCluster.push_back(ROOT::Math::XYZVector(centerOfGravity[0], centerOfGravity[1], centerOfGravity[2]));
      // Readout of the cluster weights
      for (const cell_index& cellIdx : cluster.getCells()) {
        unsigned short element = houghSpace[cellIdx[0]][cellIdx[1]][cellIdx[2]];
        readoutCluster.push_back(ROOT::Math::XYZVector(element, 0, 0));
      }
      // Readout of the cluster cell indices
      for (const cell_index& cellIdx : cluster.getCells()) {
        readoutCluster.push_back(ROOT::Math::XYZVector(cellIdx[0], cellIdx[1], cellIdx[2]));
      }
    }
    m_ndFinderTracks.push_back(NDFinderTrack(estimatedParameters,
                                             std::move(cluster), std::move(readoutHoughSpace), std::move(readoutCluster)));
  }
}

◆ transformTrackParameters()

std::array< double, 3 > transformTrackParameters ( const std::array< double, 3 > & estimatedParameters )

protected

Transform to physical units.

Definition at line 476 of file NDFinder.cc.

{
  std::array<double, 3> transformed;
  // Omega
  if (estimatedParameters[0] == 0.) {
    transformed[0] = estimatedParameters[0];
  } else { // omega = sign(q)/r, r in cm
    transformed[0] = -1 / getTrackRadius(1. / estimatedParameters[0]);
  }
  // Phi
  if (estimatedParameters[1] > 180) {
    transformed[1] = (estimatedParameters[1] - 360) * TMath::DegToRad();
  } else {
    transformed[1] = (estimatedParameters[1]) * TMath::DegToRad();
  }
  // Cot
  transformed[2] = estimatedParameters[2];
  return transformed;
}

◆ writeHitToHoughSpace()

void writeHitToHoughSpace	(	const WireInfo &	hitInfo,
		const c5array &	expHitsToWeights )

protected

Write (add) a single hit (Hough curve) to the Hough space.

Definition at line 232 of file NDFinder.cc.

{
  c3array& houghSpace = *m_houghSpace;
  for (unsigned short cotIdx = 0; cotIdx < m_nCot; ++cotIdx) {
    for (unsigned short omegaIdx = 0; omegaIdx < m_nOmega; ++omegaIdx) {
      for (unsigned short phiIdx = 0; phiIdx < m_nPhi; ++phiIdx) {
        unsigned short houghPhiIdx = (phiIdx + wireInfo.phiSectorStart) % m_nPhi;
        houghSpace[omegaIdx][houghPhiIdx][cotIdx] +=
          expHitsToWeights[wireInfo.relativeWireID][wireInfo.priorityWire][omegaIdx][phiIdx][cotIdx];
      }
    }
  }
}

Member Data Documentation

◆ m_acceptanceRanges

std::array<std::array<double, 2>, 3> m_acceptanceRanges = {m_omegaRange, m_phiRange, m_cotRange}

staticconstexprprivate

Definition at line 244 of file NDFinder.h.

244{m_omegaRange, m_phiRange, m_cotRange};

◆ m_binSizeCot

double m_binSizeCot = (m_cotRange[1] - m_cotRange[0]) / m_nCot

staticconstexprprivate

-0.377

Definition at line 243 of file NDFinder.h.

◆ m_binSizeOmega

double m_binSizeOmega = (m_omegaRange[1] - m_omegaRange[0]) / m_nOmega

staticconstexprprivate

0.2

Definition at line 241 of file NDFinder.h.

◆ m_binSizePhi

double m_binSizePhi = (m_phiRange[1] - m_phiRange[0]) / m_nPhiSector

staticconstexprprivate

0.9375

Definition at line 242 of file NDFinder.h.

◆ m_binSizes

std::array<double, 3> m_binSizes = {m_binSizeOmega, m_binSizePhi, m_binSizeCot}

staticconstexprprivate

Definition at line 245 of file NDFinder.h.

245{m_binSizeOmega, m_binSizePhi, m_binSizeCot};

◆ m_clusterer

Clusterizend m_clusterer

private

Clustering module.

Definition at line 204 of file NDFinder.h.

◆ m_compAxialBins

SectorBinning m_compAxialBins = {m_nOmega, m_nPhiComp, 1, m_nAxial, m_nPrio}

staticconstexprprivate

40, 15, 1, 41, 3

Definition at line 233 of file NDFinder.h.

233{m_nOmega, m_nPhiComp, 1, m_nAxial, m_nPrio};

◆ m_compAxialHitReps

c5array* m_compAxialHitReps = nullptr

private

m_compAxialHitReps/m_compStereoHitReps (~ Compressed in phi (width, start, values)) 5D array mapping:

[hitID]: Relative hit number of the track segment (axial [0, 40], stereo [0, 31])
[priorityWire]: Hit priority wire ([0, 2])
[omegaIdx]: Omega index of the Hough space ([0, m_nOmega - 1])
[phiIdx]: Phi start value, number of phi bins, phi values (0, 1, [2, 14])
[cotIdx]: Cot index of the Hough space (0 for axial TS, [0, m_nCot - 1] for stereo TS)

to the Hough space weight contribution at the corresponding bin (int, [0, 7])

Definition at line 267 of file NDFinder.h.

◆ m_compStereoBins

SectorBinning m_compStereoBins = {m_nOmega, m_nPhiComp, m_nCot, m_nStereo, m_nPrio}

staticconstexprprivate

40, 15, 9, 32, 3

Definition at line 234 of file NDFinder.h.

234{m_nOmega, m_nPhiComp, m_nCot, m_nStereo, m_nPrio};

◆ m_compStereoHitReps

c5array* m_compStereoHitReps = nullptr

private

Definition at line 268 of file NDFinder.h.

◆ m_cotRange

std::array<double, 2> m_cotRange = {2.3849627654510415, -1.0061730449796316}

staticconstexprprivate

=> theta in [22.75, 135.18]

Definition at line 240 of file NDFinder.h.

240{2.3849627654510415, -1.0061730449796316};

◆ m_expAxialHitReps

c5array* m_expAxialHitReps = nullptr

private

m_expAxialHitReps/m_expStereoHitReps (~ expansion of the compressed representations) 5D array mapping:

[hitID]: Relative hit number of the track segment (axial [0, 40], stereo [0, 31])
[priorityWire]: Hit priority wire ([0, 2])
[omegaIdx]: Omega index of the Hough space ([0, m_nOmega - 1])
[phiIdx]: Phi index of the Hough space ([0, m_nPhi - 1])
[cotIdx]: Cot index of the Hough space ([0, m_nCot - 1] for both axial and stereo!)

to the Hough space weight contribution at the corresponding bin (int, [0, 7])

Definition at line 280 of file NDFinder.h.

◆ m_expStereoHitReps

c5array* m_expStereoHitReps = nullptr

private

Definition at line 281 of file NDFinder.h.

◆ m_hitIDs

std::vector<unsigned short> m_hitIDs

private

TS-IDs of the hits in the current event: Elements = [0,2335] for 2336 TS in total.

Definition at line 192 of file NDFinder.h.

◆ m_hitSLIDs

std::vector<unsigned short> m_hitSLIDs

private

SL-IDs of the hits in the current event: Elements = Super layer number in [0,1,...,8].

Definition at line 194 of file NDFinder.h.

◆ m_hitToSectorIDs

c2array* m_hitToSectorIDs = nullptr

private

m_hitToSectorIDs: 2D array mapping TS-ID ([0, 2335]) to:

[0]: Orientation (1 = axial, 0 = stereo)
[1]: Relative wire ID in the sector ([0, 40] for axials, [0, 31] for stereos)
[2]: Relative phi-sector ID in the super layer ([0, 31] in each SL)

Definition at line 255 of file NDFinder.h.

◆ m_houghSpace

c3array* m_houghSpace = nullptr

private

The complete Hough space with the size [m_nOmega, m_nPhi, m_nCot].

Definition at line 283 of file NDFinder.h.

◆ m_nAxial

unsigned short m_nAxial = 41

staticconstexprprivate

Number of unique axial track segments.

Definition at line 216 of file NDFinder.h.

◆ m_nCot

unsigned short m_nCot = 9

staticconstexprprivate

Bins in the cot dimension.

Definition at line 223 of file NDFinder.h.

◆ m_ndFinderParams

NDFinderParameters m_ndFinderParams

private

Configuration parameters of the 3DFinder.

Definition at line 202 of file NDFinder.h.

◆ m_ndFinderTracks

std::vector<NDFinderTrack> m_ndFinderTracks

private

Result: Vector of the found tracks.

Definition at line 190 of file NDFinder.h.

◆ m_nHits

unsigned short m_nHits {0}

private

Counter for the number of hits in the current event.

Definition at line 200 of file NDFinder.h.

200{0};

◆ m_nOmega

unsigned short m_nOmega = 40

staticconstexprprivate

Bins in the phi dimension.

Definition at line 221 of file NDFinder.h.

◆ m_nPhi

unsigned short m_nPhi = 384

staticconstexprprivate

Bins in the omega dimension.

Definition at line 222 of file NDFinder.h.

◆ m_nPhiComp

unsigned short m_nPhiComp = 15

staticconstexprprivate

Bins of compressed phi: phi_start, phi_width, phi_0, ..., phi_12.

Definition at line 232 of file NDFinder.h.

◆ m_nPhiSector

unsigned short m_nPhiSector = m_nPhi / m_phiGeo

staticconstexprprivate

Bins of one phi sector (12)

Definition at line 229 of file NDFinder.h.

◆ m_nPrio

unsigned short m_nPrio = 3

staticconstexprprivate

Number of priority wires.

Definition at line 218 of file NDFinder.h.

◆ m_nSL

unsigned short m_nSL = 9

staticconstexprprivate

Number of super layers.

Definition at line 215 of file NDFinder.h.

◆ m_nStereo

unsigned short m_nStereo = 32

staticconstexprprivate

Number of unique stereo track segments.

Definition at line 217 of file NDFinder.h.

◆ m_nTS

unsigned short m_nTS = 2336

staticconstexprprivate

Number of track segments.

Definition at line 214 of file NDFinder.h.

◆ m_omegaRange

std::array<double, 2> m_omegaRange = {-4., 4.}

staticconstexprprivate

1/4 = 0.25 GeV (minimum transverse momentum)

Definition at line 237 of file NDFinder.h.

237{-4., 4.};

◆ m_phiGeo

unsigned short m_phiGeo = 32

staticconstexprprivate

Repetition of the wire pattern.

Definition at line 226 of file NDFinder.h.

◆ m_phiRange

std::array<double, 2> m_phiRange = {0., 11.25}

staticconstexprprivate

One phi sector (360/32)

Definition at line 238 of file NDFinder.h.

238{0., 11.25};

◆ m_priorityWirePos

std::vector<unsigned short> m_priorityWirePos

private

Priority positon within the TS. Elements basf2: [0,3] first, left, right, no hit.

Definition at line 196 of file NDFinder.h.

◆ m_priorityWireTime

std::vector<short> m_priorityWireTime

private

Drift time of the priority wire.

Definition at line 198 of file NDFinder.h.

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

trg/cdc/include/NDFinder.h
trg/cdc/src/NDFinder.cc

Classes

Public Member Functions

Protected Member Functions

Static Protected Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ ~NDFinder()

Member Function Documentation

◆ addHit()

◆ calculateCenterOfGravity()

◆ checkHitSuperLayers()

◆ extractContributionInfos()

◆ fillExpandedHitReps()

◆ findTracks()

◆ getFinderTracks()

◆ getHitContribution()

◆ getHitsVsClustersTable()

◆ getMaximumHitInSuperLayer()

◆ getTrackParameterEstimate()

◆ getTrackRadius()

◆ init()

◆ initHitToSectorMap()

◆ initLookUpArrays()

◆ loadCompressedHitReps()

◆ processHitForHoughSpace()

◆ relateHitsToClusters()

◆ reset()

◆ runTrackFinding()

◆ transformTrackParameters()

◆ writeHitToHoughSpace()

Member Data Documentation

◆ m_acceptanceRanges

◆ m_binSizeCot

◆ m_binSizeOmega

◆ m_binSizePhi

◆ m_binSizes

◆ m_clusterer

◆ m_compAxialBins

◆ m_compAxialHitReps

◆ m_compStereoBins

◆ m_compStereoHitReps

◆ m_cotRange

◆ m_expAxialHitReps

◆ m_expStereoHitReps

◆ m_hitIDs

◆ m_hitSLIDs

◆ m_hitToSectorIDs

◆ m_houghSpace

◆ m_nAxial

◆ m_nCot

◆ m_ndFinderParams

◆ m_ndFinderTracks

◆ m_nHits

◆ m_nOmega

◆ m_nPhi

◆ m_nPhiComp

◆ m_nPhiSector

◆ m_nPrio

◆ m_nSL

◆ m_nStereo

◆ m_nTS

◆ m_omegaRange

◆ m_phiGeo

◆ m_phiRange

◆ m_priorityWirePos

◆ m_priorityWireTime