LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TwoViewTransverseTracksAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
15 
17 
18 using namespace pandora;
19 
20 namespace lar_content
21 {
22 
23 TwoViewTransverseTracksAlgorithm::TwoViewTransverseTracksAlgorithm() :
24  m_nMaxMatrixToolRepeats(1000),
25  m_downsampleFactor(5),
26  m_minSamples(11),
27  m_nPermutations(1000),
28  m_localMatchingScoreThreshold(0.99f),
29  m_maxDotProduct(0.998f),
30  m_minOverallMatchingScore(0.1f),
31  m_minOverallLocallyMatchedFraction(0.1f),
32  m_randomNumberGenerator(static_cast<std::mt19937::result_type>(0))
33 {
34 }
35 
36 //------------------------------------------------------------------------------------------------------------------------------------------
37 
38 void TwoViewTransverseTracksAlgorithm::CalculateOverlapResult(const Cluster *const pCluster1, const Cluster *const pCluster2, const Cluster *const)
39 {
41  static_cast<std::mt19937::result_type>(pCluster1->GetOrderedCaloHitList().size() + pCluster2->GetOrderedCaloHitList().size()));
42 
43  TwoViewTransverseOverlapResult overlapResult;
44  PANDORA_THROW_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, this->CalculateOverlapResult(pCluster1, pCluster2, overlapResult));
45 
46  if (overlapResult.IsInitialized())
47  this->GetMatchingControl().GetOverlapMatrix().SetOverlapResult(pCluster1, pCluster2, overlapResult);
48 }
49 
50 //------------------------------------------------------------------------------------------------------------------------------------------
51 
53  const Cluster *const pCluster1, const Cluster *const pCluster2, TwoViewTransverseOverlapResult &overlapResult)
54 {
55  UIntSet daughterVolumeIntersection;
56  LArGeometryHelper::GetCommonDaughterVolumes(pCluster1, pCluster2, daughterVolumeIntersection);
57 
58  if (daughterVolumeIntersection.empty())
59  return STATUS_CODE_NOT_FOUND;
60 
61  if (this->GetPrimaryAxisDotDriftAxis(pCluster1) > m_maxDotProduct || this->GetPrimaryAxisDotDriftAxis(pCluster2) > m_maxDotProduct)
62  return STATUS_CODE_NOT_FOUND;
63 
64  float xMin1(0.f), xMax1(0.f), xMin2(0.f), xMax2(0.f);
65  pCluster1->GetClusterSpanX(xMin1, xMax1);
66  pCluster2->GetClusterSpanX(xMin2, xMax2);
67 
68  const TwoViewXOverlap twoViewXOverlap(xMin1, xMax1, xMin2, xMax2);
69  if (twoViewXOverlap.GetXSpan0() < std::numeric_limits<float>::epsilon() || twoViewXOverlap.GetXSpan1() < std::numeric_limits<float>::epsilon())
70  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
71 
72  if (twoViewXOverlap.GetTwoViewXOverlapSpan() < std::numeric_limits<float>::epsilon())
73  return STATUS_CODE_NOT_FOUND;
74 
75  float xOverlapMin(twoViewXOverlap.GetTwoViewXOverlapMin());
76  float xOverlapMax(twoViewXOverlap.GetTwoViewXOverlapMax());
77  float zMin1(0.f), zMax1(0.f);
78  float zMin2(0.f), zMax2(0.f);
79  pCluster1->GetClusterSpanZ(xMin1, xMax1, zMin1, zMax1);
80  pCluster2->GetClusterSpanZ(xMin2, xMax2, zMin2, zMax2);
81  const CartesianVector boundingBoxMin1(xOverlapMin, 0.f, zMin1), boundingBoxMax1(xOverlapMax, 0.f, zMax1);
82  const CartesianVector boundingBoxMin2(xOverlapMin, 0.f, zMin2), boundingBoxMax2(xOverlapMax, 0.f, zMax2);
83 
84  pandora::CaloHitList overlapHits1, overlapHits2;
85  LArClusterHelper::GetCaloHitListInBoundingBox(pCluster1, boundingBoxMin1, boundingBoxMax1, overlapHits1);
86  LArClusterHelper::GetCaloHitListInBoundingBox(pCluster2, boundingBoxMin2, boundingBoxMax2, overlapHits2);
87 
88  if (m_minSamples > std::min(overlapHits1.size(), overlapHits2.size()))
89  return STATUS_CODE_NOT_FOUND;
90 
91  if (1 > m_downsampleFactor)
92  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
93  const unsigned int nSamples(
94  std::max(m_minSamples, static_cast<unsigned int>(std::min(overlapHits1.size(), overlapHits2.size())) / m_downsampleFactor));
95 
97  for (const pandora::CaloHit *const pCaloHit : overlapHits1)
98  inputData1.emplace_back(pCaloHit->GetPositionVector().GetX(), pCaloHit->GetInputEnergy());
99 
101  for (const pandora::CaloHit *const pCaloHit : overlapHits2)
102  inputData2.emplace_back(pCaloHit->GetPositionVector().GetX(), pCaloHit->GetInputEnergy());
103 
104  const DiscreteProbabilityVector discreteProbabilityVector1(inputData1, xOverlapMax, false);
105  const DiscreteProbabilityVector discreteProbabilityVector2(inputData2, xOverlapMax, false);
106 
108  for (unsigned int iSample = 0; iSample < nSamples; ++iSample)
109  {
110  resamplingPointsX.emplace_back(
111  (xOverlapMin + (xOverlapMax - xOverlapMin) * static_cast<float>(iSample + 1) / static_cast<float>(nSamples + 1)));
112  }
113 
114  const DiscreteProbabilityVector resampledDiscreteProbabilityVector1(discreteProbabilityVector1, resamplingPointsX);
115  const DiscreteProbabilityVector resampledDiscreteProbabilityVector2(discreteProbabilityVector2, resamplingPointsX);
116 
117  const float correlation(
118  LArDiscreteProbabilityHelper::CalculateCorrelationCoefficient(resampledDiscreteProbabilityVector1, resampledDiscreteProbabilityVector2));
119 
121  resampledDiscreteProbabilityVector1, resampledDiscreteProbabilityVector2, m_randomNumberGenerator, m_nPermutations));
122 
123  const float matchingScore(1.f - pvalue);
124  if (matchingScore < m_minOverallMatchingScore)
125  return STATUS_CODE_NOT_FOUND;
126 
127  const unsigned int nLocallyMatchedSamplingPoints(this->CalculateNumberOfLocallyMatchingSamplingPoints(
128  resampledDiscreteProbabilityVector1, resampledDiscreteProbabilityVector2, m_randomNumberGenerator));
129  const int nComparisons(static_cast<int>(resampledDiscreteProbabilityVector1.GetSize()) - (static_cast<int>(m_minSamples) - 1));
130  if (1 > nComparisons)
131  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
132 
133  const float locallyMatchedFraction(static_cast<float>(nLocallyMatchedSamplingPoints) / static_cast<float>(nComparisons));
134  if (locallyMatchedFraction < m_minOverallLocallyMatchedFraction)
135  return STATUS_CODE_NOT_FOUND;
136 
137  overlapResult = TwoViewTransverseOverlapResult(
138  matchingScore, m_downsampleFactor, nComparisons, nLocallyMatchedSamplingPoints, correlation, twoViewXOverlap);
139 
140  return STATUS_CODE_SUCCESS;
141 }
142 
143 //------------------------------------------------------------------------------------------------------------------------------------------
144 
146  const DiscreteProbabilityVector &discreteProbabilityVector2, std::mt19937 &randomNumberGenerator)
147 {
148  if (discreteProbabilityVector1.GetSize() != discreteProbabilityVector2.GetSize() ||
149  0 == discreteProbabilityVector1.GetSize() * discreteProbabilityVector2.GetSize())
150  throw STATUS_CODE_INVALID_PARAMETER;
151 
152  if (m_minSamples > discreteProbabilityVector1.GetSize())
153  throw STATUS_CODE_INVALID_PARAMETER;
154 
155  pandora::FloatVector localValues1, localValues2;
156  unsigned int nMatchedComparisons(0);
157 
158  for (unsigned int iValue = 0; iValue < discreteProbabilityVector1.GetSize(); ++iValue)
159  {
160  localValues1.emplace_back(discreteProbabilityVector1.GetProbability(iValue));
161  localValues2.emplace_back(discreteProbabilityVector2.GetProbability(iValue));
162  if (localValues1.size() == m_minSamples)
163  {
164  float localPValue(0);
165  try
166  {
168  localValues1, localValues2, randomNumberGenerator, m_nPermutations);
169  }
170  catch (const StatusCodeException &)
171  {
172  std::cout << "TwoViewTransverseTracksAlgorithm: failed to calculate correlation coefficient p-value for these numbers" << std::endl;
173  ;
174  std::cout << "----view 0: ";
175  for (unsigned int iElement = 0; iElement < localValues1.size(); ++iElement)
176  std::cout << localValues1.at(iElement) << " ";
177  std::cout << std::endl;
178  std::cout << "----view 1: ";
179  for (unsigned int iElement = 0; iElement < localValues2.size(); ++iElement)
180  std::cout << localValues2.at(iElement) << " ";
181  std::cout << std::endl;
182  }
183 
184  if ((1.f - localPValue) - m_localMatchingScoreThreshold > std::numeric_limits<float>::epsilon())
185  nMatchedComparisons++;
186 
187  localValues1.erase(localValues1.begin());
188  localValues2.erase(localValues2.begin());
189  }
190  }
191  return nMatchedComparisons;
192 }
193 
194 //------------------------------------------------------------------------------------------------------------------------------------------
195 
196 float TwoViewTransverseTracksAlgorithm::GetPrimaryAxisDotDriftAxis(const pandora::Cluster *const pCluster)
197 {
198  pandora::CartesianPointVector pointVector;
199  LArClusterHelper::GetCoordinateVector(pCluster, pointVector);
200 
201  pandora::CartesianVector centroid(0.f, 0.f, 0.f);
202  LArPcaHelper::EigenVectors eigenVecs;
203  LArPcaHelper::EigenValues eigenValues(0.f, 0.f, 0.f);
204  LArPcaHelper::RunPca(pointVector, centroid, eigenValues, eigenVecs);
205 
206  const pandora::CartesianVector primaryAxis(eigenVecs.at(0));
207  const pandora::CartesianVector driftAxis(1.f, 0.f, 0.f);
208  return primaryAxis.GetDotProduct(driftAxis);
209 }
210 
211 //------------------------------------------------------------------------------------------------------------------------------------------
212 
214 {
215  unsigned int repeatCounter(0);
216  for (MatrixToolVector::const_iterator iter = m_algorithmToolVector.begin(), iterEnd = m_algorithmToolVector.end(); iter != iterEnd;)
217  {
218  if ((*iter)->Run(this, this->GetMatchingControl().GetOverlapMatrix()))
219  {
220  iter = m_algorithmToolVector.begin();
221 
222  if (++repeatCounter > m_nMaxMatrixToolRepeats)
223  break;
224  }
225  else
226  {
227  ++iter;
228  }
229  }
230 }
231 
232 //------------------------------------------------------------------------------------------------------------------------------------------
233 
234 StatusCode TwoViewTransverseTracksAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
235 {
236  AlgorithmToolVector algorithmToolVector;
237  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ProcessAlgorithmToolList(*this, xmlHandle, "TrackTools", algorithmToolVector));
238 
239  for (AlgorithmToolVector::const_iterator iter = algorithmToolVector.begin(), iterEnd = algorithmToolVector.end(); iter != iterEnd; ++iter)
240  {
241  TransverseMatrixTool *const pTransverseMatrixTool(dynamic_cast<TransverseMatrixTool *>(*iter));
242 
243  if (!pTransverseMatrixTool)
244  return STATUS_CODE_INVALID_PARAMETER;
245 
246  m_algorithmToolVector.push_back(pTransverseMatrixTool);
247  }
248 
249  PANDORA_RETURN_RESULT_IF_AND_IF(
250  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "NMaxMatrixToolRepeats", m_nMaxMatrixToolRepeats));
251 
252  PANDORA_RETURN_RESULT_IF_AND_IF(
253  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DownsampleFactor", m_downsampleFactor));
254 
255  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinSamples", m_minSamples));
256 
257  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "NPermutations", m_nPermutations));
258 
259  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
260  XmlHelper::ReadValue(xmlHandle, "LocalMatchingScoreThreshold", m_localMatchingScoreThreshold));
261 
262  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxDotProduct", m_maxDotProduct));
263 
264  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
265  XmlHelper::ReadValue(xmlHandle, "MinOverallMatchingScore", m_minOverallMatchingScore));
266 
267  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
268  XmlHelper::ReadValue(xmlHandle, "MinOverallLocallyMatchedFraction", m_minOverallLocallyMatchedFraction));
269 
270  return BaseAlgorithm::ReadSettings(xmlHandle);
271 }
272 
273 } // namespace lar_content
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
bool IsInitialized() const
Whether the track overlap result has been initialized.
std::mt19937 m_randomNumberGenerator
The random number generator.
float GetTwoViewXOverlapMax() const
Get the x overlap min X value.
pandora::CartesianVector EigenValues
Definition: LArPcaHelper.h:24
MatrixType & GetOverlapMatrix()
Get the overlap matrix.
STL namespace.
static float CalculateCorrelationCoefficientPValueFromPermutationTest(const T &t1, const T &t2, std::mt19937 &randomNumberGenerator, const unsigned int nPermutations)
Calculate P value for measured correlation coefficient between two datasets via a permutation test...
Header file for the principal curve analysis helper class.
intermediate_table::const_iterator const_iterator
float m_localMatchingScoreThreshold
The minimum score to classify a local region as matching.
MatrixToolVector m_algorithmToolVector
The algorithm tool vector.
static float CalculateCorrelationCoefficient(const T &t1, const T &t2)
Calculate the correlation coefficient between two datasets.
Header file for the discrete probability helper class.
TFile f
Definition: plotHisto.C:6
float GetXSpan0() const
Get the x span in the view 0.
Header file for the geometry helper class.
float GetTwoViewXOverlapMin() const
Get the x overlap max X value.
float GetXSpan1() const
Get the x span in the view 1.
Header file for the cluster helper class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
static void GetCaloHitListInBoundingBox(const pandora::Cluster *const pCluster, const pandora::CartesianVector &lowerBound, const pandora::CartesianVector &upperBound, pandora::CaloHitList &caloHitList)
Get list of Calo hits from an input cluster that are contained in a bounding box. The hits are sorted...
float GetProbability(const unsigned int index) const
Get the probability value of the element in the vector.
float GetPrimaryAxisDotDriftAxis(const pandora::Cluster *const pCluster)
Get the dot product between the cluster&#39;s primary axis and the drift axis.
unsigned int m_nPermutations
The number of permutations for calculating p-values.
Header file for the two view transverse tracks algorithm class.
void SetOverlapResult(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2, const OverlapResult &overlapResult)
Set overlap result.
unsigned int CalculateNumberOfLocallyMatchingSamplingPoints(const DiscreteProbabilityVector &discreteProbabilityVector1, const DiscreteProbabilityVector &discreteProbabilityVector2, std::mt19937 &randomNumberGenerator)
Calculates the number of the sliding windows that contains charge bins that locally match...
unsigned int m_minSamples
The minimum number of samples needed for comparing charges.
float m_minOverallMatchingScore
M The maximum allowed cluster primary qxis Dot drift axis to fill the overlap result.
float m_minOverallLocallyMatchedFraction
The minimum required lcoally matched fraction to fill the overlap result.
static void RunPca(const T &t, pandora::CartesianVector &centroid, EigenValues &outputEigenValues, EigenVectors &outputEigenVectors)
Run principal component analysis using input calo hits (TPC_VIEW_U,V,W or TPC_3D; all treated as 3D p...
void ExamineOverlapContainer()
Examine contents of overlap container, collect together best-matching 2D particles and modify cluster...
std::vector< pandora::CartesianVector > EigenVectors
Definition: LArPcaHelper.h:25
unsigned int m_downsampleFactor
The downsampling (hit merging) applied to hits in the overlap region.
float GetTwoViewXOverlapSpan() const
Get the x overlap span.
unsigned int GetSize() const
Get the size of the probability vector.
static void GetCoordinateVector(const pandora::Cluster *const pCluster, pandora::CartesianPointVector &coordinateVector)
Get vector of hit coordinates from an input cluster.
TwoViewXOverlap class.
void CalculateOverlapResult(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2, const pandora::Cluster *const)
Calculate cluster overlap result and store in container.
unsigned int m_nMaxMatrixToolRepeats
The maximum number of repeat loops over matrix tools.
static void GetCommonDaughterVolumes(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2, UIntSet &intersect)
Return the set of common daughter volumes between two 2D clusters.