Pandora
Pandora source code navigator
Loading...
Searching...
No Matches
ParticleRecoveryAlgorithm.cc
Go to the documentation of this file.
1
10
14
16
18
19#include <algorithm>
20
21using namespace pandora;
22
23namespace lar_content
24{
25
27 m_checkGaps(true),
28 m_minClusterCaloHits(20),
29 m_minClusterLengthSquared(5.f * 5.f),
30 m_minClusterXSpan(0.25f),
31 m_vertexClusterMode(false),
32 m_minVertexLongitudinalDistance(-2.5f),
33 m_maxVertexTransverseDistance(1.5f),
34 m_minXOverlapFraction(0.75f),
35 m_minXOverlapFractionGaps(0.75f),
36 m_sampleStepSize(0.5f),
37 m_slidingFitHalfWindow(10),
38 m_pseudoChi2Cut(5.f)
39{
40}
41
42//------------------------------------------------------------------------------------------------------------------------------------------
43
45{
46 ClusterList inputClusterListU, inputClusterListV, inputClusterListW;
47 this->GetInputClusters(inputClusterListU, inputClusterListV, inputClusterListW);
48
49 ClusterList selectedClusterListU, selectedClusterListV, selectedClusterListW;
50 this->SelectInputClusters(inputClusterListU, selectedClusterListU);
51 this->SelectInputClusters(inputClusterListV, selectedClusterListV);
52 this->SelectInputClusters(inputClusterListW, selectedClusterListW);
53
54 SimpleOverlapTensor overlapTensor;
55 this->FindOverlaps(selectedClusterListU, selectedClusterListV, overlapTensor);
56 this->FindOverlaps(selectedClusterListV, selectedClusterListW, overlapTensor);
57 this->FindOverlaps(selectedClusterListW, selectedClusterListU, overlapTensor);
58 this->ExamineTensor(overlapTensor);
59
60 return STATUS_CODE_SUCCESS;
61}
62
63//------------------------------------------------------------------------------------------------------------------------------------------
64
65void ParticleRecoveryAlgorithm::GetInputClusters(ClusterList &inputClusterListU, ClusterList &inputClusterListV, ClusterList &inputClusterListW) const
66{
67 for (StringVector::const_iterator iter = m_inputClusterListNames.begin(), iterEnd = m_inputClusterListNames.end(); iter != iterEnd; ++iter)
68 {
69 const ClusterList *pClusterList(NULL);
70 PANDORA_THROW_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, *iter, pClusterList));
71
72 if (!pClusterList || pClusterList->empty())
73 {
75 std::cout << "ParticleRecoveryAlgorithm: unable to find cluster list " << *iter << std::endl;
76
77 continue;
78 }
79
80 for (ClusterList::const_iterator cIter = pClusterList->begin(), cIterEnd = pClusterList->end(); cIter != cIterEnd; ++cIter)
81 {
82 const Cluster *const pCluster(*cIter);
83 const HitType hitType(LArClusterHelper::GetClusterHitType(pCluster));
84
85 if ((TPC_VIEW_U != hitType) && (TPC_VIEW_V != hitType) && (TPC_VIEW_W != hitType))
86 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
87
88 ClusterList &clusterList((TPC_VIEW_U == hitType) ? inputClusterListU : (TPC_VIEW_V == hitType) ? inputClusterListV : inputClusterListW);
89 clusterList.push_back(pCluster);
90 }
91 }
92}
93
94//------------------------------------------------------------------------------------------------------------------------------------------
95
96void ParticleRecoveryAlgorithm::SelectInputClusters(const ClusterList &inputClusterList, ClusterList &selectedClusterList) const
97{
99 {
100 ClusterList vertexClusterList;
101 this->VertexClusterSelection(inputClusterList, vertexClusterList);
102 this->StandardClusterSelection(vertexClusterList, selectedClusterList);
103 }
104 else
105 {
106 this->StandardClusterSelection(inputClusterList, selectedClusterList);
107 }
108}
109
110//------------------------------------------------------------------------------------------------------------------------------------------
111
112void ParticleRecoveryAlgorithm::StandardClusterSelection(const ClusterList &inputClusterList, ClusterList &selectedClusterList) const
113{
114 for (ClusterList::const_iterator iter = inputClusterList.begin(), iterEnd = inputClusterList.end(); iter != iterEnd; ++iter)
115 {
116 const Cluster *const pCluster = *iter;
117
118 if (!pCluster->IsAvailable())
119 continue;
120
121 if (pCluster->GetNCaloHits() < m_minClusterCaloHits)
122 continue;
123
125 continue;
126
127 float xMin(0.f), xMax(0.f);
128 pCluster->GetClusterSpanX(xMin, xMax);
129
130 if ((xMax - xMin) < m_minClusterXSpan)
131 continue;
132
133 selectedClusterList.push_back(pCluster);
134 }
135}
136
137//------------------------------------------------------------------------------------------------------------------------------------------
138
139void ParticleRecoveryAlgorithm::VertexClusterSelection(const ClusterList &inputClusterList, ClusterList &selectedClusterList) const
140{
141 CartesianPointVector vertexList;
142
143 for (ClusterList::const_iterator iter = inputClusterList.begin(), iterEnd = inputClusterList.end(); iter != iterEnd; ++iter)
144 {
145 try
146 {
147 if (!(*iter)->IsAvailable())
148 {
149 const LArPointingCluster pointingCluster(*iter);
150 vertexList.push_back(pointingCluster.GetInnerVertex().GetPosition());
151 vertexList.push_back(pointingCluster.GetOuterVertex().GetPosition());
152 }
153 }
154 catch (StatusCodeException &)
155 {
156 }
157 }
158
159 for (ClusterList::const_iterator iter = inputClusterList.begin(), iterEnd = inputClusterList.end(); iter != iterEnd; ++iter)
160 {
161 try
162 {
163 const Cluster *const pCluster = *iter;
164
165 if (!pCluster->IsAvailable())
166 continue;
167
168 const LArPointingCluster pointingCluster(pCluster);
169
170 for (CartesianPointVector::const_iterator vIter = vertexList.begin(), vIterEnd = vertexList.end(); vIter != vIterEnd; ++vIter)
171 {
174 {
175 selectedClusterList.push_back(pCluster);
176 break;
177 }
178 }
179 }
180 catch (StatusCodeException &)
181 {
182 }
183 }
184}
185
186//------------------------------------------------------------------------------------------------------------------------------------------
187
188void ParticleRecoveryAlgorithm::FindOverlaps(const ClusterList &clusterList1, const ClusterList &clusterList2, SimpleOverlapTensor &overlapTensor) const
189{
190 for (ClusterList::const_iterator iter1 = clusterList1.begin(), iter1End = clusterList1.end(); iter1 != iter1End; ++iter1)
191 {
192 for (ClusterList::const_iterator iter2 = clusterList2.begin(), iter2End = clusterList2.end(); iter2 != iter2End; ++iter2)
193 {
194 if (this->IsOverlap(*iter1, *iter2))
195 overlapTensor.AddAssociation(*iter1, *iter2);
196 }
197 }
198}
199
200//------------------------------------------------------------------------------------------------------------------------------------------
201
202bool ParticleRecoveryAlgorithm::IsOverlap(const Cluster *const pCluster1, const Cluster *const pCluster2) const
203{
205 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
206
207 if ((0 == pCluster1->GetNCaloHits()) || (0 == pCluster2->GetNCaloHits()))
208 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
209
210 float xMin1(0.f), xMax1(0.f), xMin2(0.f), xMax2(0.f);
211 pCluster1->GetClusterSpanX(xMin1, xMax1);
212 pCluster2->GetClusterSpanX(xMin2, xMax2);
213
214 const float xSpan1(xMax1 - xMin1), xSpan2(xMax2 - xMin2);
215
216 if ((xSpan1 < std::numeric_limits<float>::epsilon()) || (xSpan2 < std::numeric_limits<float>::epsilon()))
217 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
218
219 const float xOverlap(std::min(xMax1, xMax2) - std::max(xMin1, xMin2));
220
221 float xOverlapFraction1(xOverlap / xSpan1), xOverlapFraction2(xOverlap / xSpan2);
222
223 if (m_checkGaps)
224 this->CalculateEffectiveOverlapFractions(pCluster1, xMin1, xMax1, pCluster2, xMin2, xMax2, xOverlapFraction1, xOverlapFraction2);
225
226 if ((xOverlapFraction1 < m_minXOverlapFraction) || (xOverlapFraction2 < m_minXOverlapFraction))
227 return false;
228
229 return true;
230}
231
232//------------------------------------------------------------------------------------------------------------------------------------------
233
234void ParticleRecoveryAlgorithm::CalculateEffectiveOverlapFractions(const Cluster *const pCluster1, const float xMin1, const float xMax1,
235 const Cluster *const pCluster2, const float xMin2, const float xMax2, float &xOverlapFraction1, float &xOverlapFraction2) const
236{
238 return;
239
240 const float xMin(std::min(xMin1, xMin2));
241 const float xMax(std::max(xMax1, xMax2));
242 float xMinEff1(xMin1), xMaxEff1(xMax1), xMinEff2(xMin2), xMaxEff2(xMax2);
243
244 this->CalculateEffectiveSpan(pCluster1, xMin, xMax, xMinEff1, xMaxEff1);
245 this->CalculateEffectiveSpan(pCluster2, xMin, xMax, xMinEff2, xMaxEff2);
246
247 const float effectiveXSpan1(xMaxEff1 - xMinEff1), effectiveXSpan2(xMaxEff2 - xMinEff2);
248 const float effectiveXOverlapSpan(std::min(xMaxEff1, xMaxEff2) - std::max(xMinEff1, xMinEff2));
249
250 if ((effectiveXSpan1 > std::numeric_limits<float>::epsilon()) && (effectiveXSpan2 > std::numeric_limits<float>::epsilon()))
251 {
252 xOverlapFraction1 = std::min(1.f, (effectiveXOverlapSpan / effectiveXSpan1));
253 xOverlapFraction2 = std::min(1.f, (effectiveXOverlapSpan / effectiveXSpan2));
254 }
255}
256
257//------------------------------------------------------------------------------------------------------------------------------------------
258
260 const pandora::Cluster *const pCluster, const float xMin, const float xMax, float &xMinEff, float &xMaxEff) const
261{
262 // TODO cache sliding linear fit results and optimise protection against exceptions from TwoDSlidingFitResult and IsXSamplingPointInGap
263 try
264 {
265 const float slidingFitPitch(LArGeometryHelper::GetWireZPitch(this->GetPandora()));
266
267 const TwoDSlidingFitResult slidingFitResult(pCluster, m_slidingFitHalfWindow, slidingFitPitch);
268
269 const int nSamplingPointsLeft(1 + static_cast<int>((xMinEff - xMin) / m_sampleStepSize));
270 const int nSamplingPointsRight(1 + static_cast<int>((xMax - xMaxEff) / m_sampleStepSize));
271 float dxMin(0.f), dxMax(0.f);
272
273 for (int iSample = 1; iSample <= nSamplingPointsLeft; ++iSample)
274 {
275 const float xSample(std::max(xMin, xMinEff - static_cast<float>(iSample) * m_sampleStepSize));
276
277 if (!LArGeometryHelper::IsXSamplingPointInGap(this->GetPandora(), xSample, slidingFitResult, m_sampleStepSize))
278 break;
279
280 dxMin = xMinEff - xSample;
281 }
282
283 for (int iSample = 1; iSample <= nSamplingPointsRight; ++iSample)
284 {
285 const float xSample(std::min(xMax, xMaxEff + static_cast<float>(iSample) * m_sampleStepSize));
286
287 if (!LArGeometryHelper::IsXSamplingPointInGap(this->GetPandora(), xSample, slidingFitResult, m_sampleStepSize))
288 break;
289
290 dxMax = xSample - xMaxEff;
291 }
292
293 xMinEff = xMinEff - dxMin;
294 xMaxEff = xMaxEff + dxMax;
295 }
296 catch (StatusCodeException &)
297 {
298 }
299}
300
301//------------------------------------------------------------------------------------------------------------------------------------------
302
304{
305 ClusterVector sortedKeyClusters(overlapTensor.GetKeyClusters().begin(), overlapTensor.GetKeyClusters().end());
306 std::sort(sortedKeyClusters.begin(), sortedKeyClusters.end(), LArClusterHelper::SortByNHits);
307
308 for (const Cluster *const pKeyCluster : sortedKeyClusters)
309 {
310 ClusterList clusterListU, clusterListV, clusterListW;
311
312 overlapTensor.GetConnectedElements(pKeyCluster, true, clusterListU, clusterListV, clusterListW);
313 const unsigned int nU(clusterListU.size()), nV(clusterListV.size()), nW(clusterListW.size());
314
315 if ((0 == nU * nV) && (0 == nV * nW) && (0 == nW * nU))
316 continue;
317
318 ClusterList clusterListAll;
319 clusterListAll.insert(clusterListAll.end(), clusterListU.begin(), clusterListU.end());
320 clusterListAll.insert(clusterListAll.end(), clusterListV.begin(), clusterListV.end());
321 clusterListAll.insert(clusterListAll.end(), clusterListW.begin(), clusterListW.end());
322
323 if ((1 == nU * nV * nW) && this->CheckConsistency(*(clusterListU.begin()), *(clusterListV.begin()), *(clusterListW.begin())))
324 {
325 this->CreateTrackParticle(clusterListAll);
326 }
327 else if ((0 == nU * nV * nW) && ((1 == nU && 1 == nV) || (1 == nV && 1 == nW) || (1 == nW && 1 == nU)))
328 {
329 this->CreateTrackParticle(clusterListAll);
330 }
331 else
332 {
333 // TODO - check here whether there is a gap in the 2 in one view when 1:2:0
334 // May later choose to resolve simple ambiguities, e.g. of form nU:nV:nW == 1:2:0
335 }
336 }
337}
338
339//------------------------------------------------------------------------------------------------------------------------------------------
340
341bool ParticleRecoveryAlgorithm::CheckConsistency(const Cluster *const pClusterU, const Cluster *const pClusterV, const Cluster *const pClusterW) const
342{
343 // Requirements on X matching
344 float xMinU(0.f), xMinV(0.f), xMinW(0.f), xMaxU(0.f), xMaxV(0.f), xMaxW(0.f);
345 pClusterU->GetClusterSpanX(xMinU, xMaxU);
346 pClusterV->GetClusterSpanX(xMinV, xMaxV);
347 pClusterW->GetClusterSpanX(xMinW, xMaxW);
348
349 const float xMin(std::max(xMinU, std::max(xMinV, xMinW)));
350 const float xMax(std::min(xMaxU, std::min(xMaxV, xMaxW)));
351 const float xOverlap(xMax - xMin);
352
353 if (xOverlap < std::numeric_limits<float>::epsilon())
354 return false;
355
356 // Requirements on 3D matching
357 float u(std::numeric_limits<float>::max()), v(std::numeric_limits<float>::max()), w(std::numeric_limits<float>::max());
358
359 if ((STATUS_CODE_SUCCESS != LArClusterHelper::GetAverageZ(pClusterU, xMin, xMax, u)) ||
360 (STATUS_CODE_SUCCESS != LArClusterHelper::GetAverageZ(pClusterV, xMin, xMax, v)) ||
361 (STATUS_CODE_SUCCESS != LArClusterHelper::GetAverageZ(pClusterW, xMin, xMax, w)))
362 {
363 return false;
364 }
365
366 const float uv2w(LArGeometryHelper::MergeTwoPositions(this->GetPandora(), TPC_VIEW_U, TPC_VIEW_V, u, v));
367 const float vw2u(LArGeometryHelper::MergeTwoPositions(this->GetPandora(), TPC_VIEW_V, TPC_VIEW_W, v, w));
368 const float wu2v(LArGeometryHelper::MergeTwoPositions(this->GetPandora(), TPC_VIEW_W, TPC_VIEW_U, w, u));
369
370 const float pseudoChi2(((u - vw2u) * (u - vw2u) + (v - wu2v) * (v - wu2v) + (w - uv2w) * (w - uv2w)) / 3.f);
371
372 if (pseudoChi2 > m_pseudoChi2Cut)
373 return false;
374
375 return true;
376}
377
378//------------------------------------------------------------------------------------------------------------------------------------------
379
381{
382 if (clusterList.size() < 2)
383 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
384
385 const PfoList *pPfoList = NULL;
386 std::string pfoListName;
387 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pPfoList, pfoListName));
388
389 // TODO Correct these placeholder parameters
390 PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
391 pfoParameters.m_particleId = MU_MINUS;
392 pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
393 pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
394 pfoParameters.m_energy = 0.f;
395 pfoParameters.m_momentum = CartesianVector(0.f, 0.f, 0.f);
396 pfoParameters.m_clusterList = clusterList;
397
398 const ParticleFlowObject *pPfo(NULL);
399 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pPfo));
400
401 if (!pPfoList->empty())
402 {
405 }
406}
407
408//------------------------------------------------------------------------------------------------------------------------------------------
409//------------------------------------------------------------------------------------------------------------------------------------------
410
411void ParticleRecoveryAlgorithm::SimpleOverlapTensor::AddAssociation(const Cluster *const pCluster1, const Cluster *const pCluster2)
412{
413 const HitType hitType1(LArClusterHelper::GetClusterHitType(pCluster1));
414 const HitType hitType2(LArClusterHelper::GetClusterHitType(pCluster2));
415
416 const Cluster *const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 : NULL);
417 const Cluster *const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 : NULL);
418 const Cluster *const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 : NULL);
419
420 if (pClusterU && pClusterV && !pClusterW)
421 {
422 m_clusterNavigationMapUV[pClusterU].push_back(pClusterV);
423
424 if (m_keyClusters.end() == std::find(m_keyClusters.begin(), m_keyClusters.end(), pClusterU))
425 m_keyClusters.push_back(pClusterU);
426 }
427 else if (!pClusterU && pClusterV && pClusterW)
428 {
429 m_clusterNavigationMapVW[pClusterV].push_back(pClusterW);
430
431 if (m_keyClusters.end() == std::find(m_keyClusters.begin(), m_keyClusters.end(), pClusterV))
432 m_keyClusters.push_back(pClusterV);
433 }
434 else if (pClusterU && !pClusterV && pClusterW)
435 {
436 m_clusterNavigationMapWU[pClusterW].push_back(pClusterU);
437
438 if (m_keyClusters.end() == std::find(m_keyClusters.begin(), m_keyClusters.end(), pClusterW))
439 m_keyClusters.push_back(pClusterW);
440 }
441 else
442 {
443 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
444 }
445}
446
447//------------------------------------------------------------------------------------------------------------------------------------------
448
449void ParticleRecoveryAlgorithm::SimpleOverlapTensor::GetConnectedElements(const Cluster *const pCluster, const bool ignoreUnavailable,
450 ClusterList &clusterListU, ClusterList &clusterListV, ClusterList &clusterListW) const
451{
452 if (ignoreUnavailable && !pCluster->IsAvailable())
453 return;
454
455 const HitType hitType(LArClusterHelper::GetClusterHitType(pCluster));
456
457 if (!((TPC_VIEW_U == hitType) || (TPC_VIEW_V == hitType) || (TPC_VIEW_W == hitType)))
458 throw StatusCodeException(STATUS_CODE_FAILURE);
459
460 ClusterList &clusterList((TPC_VIEW_U == hitType) ? clusterListU : (TPC_VIEW_V == hitType) ? clusterListV : clusterListW);
461 const ClusterNavigationMap &navigationMap(
462 (TPC_VIEW_U == hitType) ? m_clusterNavigationMapUV : (TPC_VIEW_V == hitType) ? m_clusterNavigationMapVW : m_clusterNavigationMapWU);
463
464 if (clusterList.end() != std::find(clusterList.begin(), clusterList.end(), pCluster))
465 return;
466
467 clusterList.push_back(pCluster);
468
469 ClusterNavigationMap::const_iterator iter = navigationMap.find(pCluster);
470
471 if (navigationMap.end() == iter)
472 return;
473
474 for (ClusterList::const_iterator cIter = iter->second.begin(), cIterEnd = iter->second.end(); cIter != cIterEnd; ++cIter)
475 this->GetConnectedElements(*cIter, ignoreUnavailable, clusterListU, clusterListV, clusterListW);
476}
477
478//------------------------------------------------------------------------------------------------------------------------------------------
479//------------------------------------------------------------------------------------------------------------------------------------------
480
482{
483 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "InputClusterListNames", m_inputClusterListNames));
484 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputPfoListName", m_outputPfoListName));
485
487 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterCaloHits", m_minClusterCaloHits));
488
489 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "CheckGaps", m_checkGaps));
490
491 float minClusterLength = std::sqrt(m_minClusterLengthSquared);
492 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterLength", minClusterLength));
493 m_minClusterLengthSquared = minClusterLength * minClusterLength;
494
495 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterXSpan", m_minClusterXSpan));
496
498 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "VertexClusterMode", m_vertexClusterMode));
499
500 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
501 XmlHelper::ReadValue(xmlHandle, "MinVertexLongitudinalDistance", m_minVertexLongitudinalDistance));
502
503 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
504 XmlHelper::ReadValue(xmlHandle, "MaxVertexTransverseDistance", m_maxVertexTransverseDistance));
505
507 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinXOverlapFraction", m_minXOverlapFraction));
508
509 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
510 XmlHelper::ReadValue(xmlHandle, "MinXOverlapFractionGaps", m_minXOverlapFractionGaps));
511
512 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SampleStepSize", m_sampleStepSize));
513
514 if (m_sampleStepSize < std::numeric_limits<float>::epsilon())
515 {
516 std::cout << "ParticleRecoveryAlgorithm: Invalid value for SampleStepSize " << m_sampleStepSize << std::endl;
517 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
518 }
519
521 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SlidingFitHalfWindow", m_slidingFitHalfWindow));
522
523 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "PseudoChi2Cut", m_pseudoChi2Cut));
524
525 return STATUS_CODE_SUCCESS;
526}
527
528} // namespace lar_content
Grouping of header files for many classes of use in particle flow algorithms.
Header file for the cluster helper class.
Header file for the geometry helper class.
Header file for the lar pointing cluster class.
Header file for the track recovery algorithm class.
#define PANDORA_THROW_RESULT_IF_AND_IF(StatusCode1, StatusCode2, Operator, Command)
Definition StatusCodes.h:55
#define PANDORA_THROW_RESULT_IF(StatusCode1, Operator, Command)
Definition StatusCodes.h:43
#define PANDORA_RETURN_RESULT_IF_AND_IF(StatusCode1, StatusCode2, Operator, Command)
Definition StatusCodes.h:31
#define PANDORA_RETURN_RESULT_IF(StatusCode1, Operator, Command)
Definition StatusCodes.h:19
static const pandora::GeometryManager * GetGeometry(const pandora::Algorithm &algorithm)
Get the pandora geometry instance.
static pandora::StatusCode ReplaceCurrentList(const pandora::Algorithm &algorithm, const std::string &newListName)
Replace the current list with a pre-saved list; use this new list as a permanent replacement for the ...
static pandora::StatusCode CreateTemporaryListAndSetCurrent(const pandora::Algorithm &algorithm, const T *&pT, std::string &temporaryListName)
Create a temporary list and set it to be the current list, enabling object creation.
static pandora::StatusCode GetList(const pandora::Algorithm &algorithm, const std::string &listName, const T *&pT)
Get a named list.
static const pandora::PandoraSettings * GetSettings(const pandora::Algorithm &algorithm)
Get the pandora settings instance.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
static pandora::StatusCode GetAverageZ(const pandora::Cluster *const pCluster, const float xmin, const float xmax, float &averageZ)
Get average Z positions of the calo hits in a cluster in range xmin to xmax.
static bool SortByNHits(const pandora::Cluster *const pLhs, const pandora::Cluster *const pRhs)
Sort clusters by number of hits, then layer span, then inner layer, then position,...
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
static float MergeTwoPositions(const pandora::Pandora &pandora, const pandora::HitType view1, const pandora::HitType view2, const float position1, const float position2)
Merge two views (U,V) to give a third view (Z).
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
static bool IsXSamplingPointInGap(const pandora::Pandora &pandora, const float xSample, const TwoDSlidingFitResult &slidingFitResult, const float gapTolerance=0.f)
Whether there is a gap in a cluster (described via its sliding fit result) at a specified x sampling ...
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
static bool IsNode(const pandora::CartesianVector &parentVertex, const LArPointingCluster::Vertex &daughterVertex, const float minLongitudinalDistance, const float maxTransverseDistance)
Whether pointing vertex is adjacent to a given position.
LArPointingCluster class.
const Vertex & GetInnerVertex() const
Get the inner vertex.
const Vertex & GetOuterVertex() const
Get the outer vertex.
ClusterNavigationMap m_clusterNavigationMapUV
The cluster navigation map U->V.
void GetConnectedElements(const pandora::Cluster *const pCluster, const bool ignoreUnavailable, pandora::ClusterList &clusterListU, pandora::ClusterList &clusterListV, pandora::ClusterList &clusterListW) const
Get elements connected to a specified cluster.
ClusterNavigationMap m_clusterNavigationMapWU
The cluster navigation map W->U.
std::unordered_map< const pandora::Cluster *, pandora::ClusterList > ClusterNavigationMap
void AddAssociation(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2)
Add an association between two clusters to the simple overlap tensor.
ClusterNavigationMap m_clusterNavigationMapVW
The cluster navigation map V->W.
const pandora::ClusterList & GetKeyClusters() const
Get the list of key clusters.
pandora::ClusterList m_keyClusters
The list of key clusters.
float m_minClusterXSpan
The min x span required in order to consider a cluster.
std::string m_outputPfoListName
The output pfo list name.
unsigned int m_minClusterCaloHits
The min number of hits in base cluster selection method.
bool m_vertexClusterMode
Whether to demand clusters are associated with vertices of existing particles.
void CreateTrackParticle(const pandora::ClusterList &clusterList) const
Create and save a track particle containing the provided clusters.
pandora::StringVector m_inputClusterListNames
The list of cluster list names.
float m_minClusterLengthSquared
The min length (squared) in base cluster selection method.
void ExamineTensor(const SimpleOverlapTensor &overlapTensor) const
Identify unambiguous cluster overlaps and resolve ambiguous overlaps, creating new track particles.
unsigned int m_slidingFitHalfWindow
The half window for the fit sliding result constructor.
void GetInputClusters(pandora::ClusterList &inputClusterListU, pandora::ClusterList &inputClusterListV, pandora::ClusterList &inputClusterListW) const
Get the input cluster lists for processing in this algorithm.
void VertexClusterSelection(const pandora::ClusterList &inputClusterList, pandora::ClusterList &selectedClusterList) const
Select a subset of input clusters nodally associated with the vertices of existing particles.
float m_minXOverlapFractionGaps
The min x overlap fraction when there are gaps involved.
void StandardClusterSelection(const pandora::ClusterList &inputClusterList, pandora::ClusterList &selectedClusterList) const
Select a subset of input clusters for processing in this algorithm.
bool m_checkGaps
Whether to check for gaps in the calculation of the overlap.
float m_maxVertexTransverseDistance
Vertex association check: max transverse distance cut.
float m_sampleStepSize
The sampling step size used in association checks, units cm.
float m_pseudoChi2Cut
The selection cut on the matched chi2.
void CalculateEffectiveOverlapFractions(const pandora::Cluster *const pCluster1, const float xMin1, const float xMax1, const pandora::Cluster *const pCluster2, const float xMin2, const float xMax2, float &xOverlapFraction1, float &xOverlapFraction2) const
Calculate effective overlap fractions taking into account gaps.
pandora::StatusCode Run()
Run the algorithm.
bool CheckConsistency(const pandora::Cluster *const pClusterU, const pandora::Cluster *const pClusterV, const pandora::Cluster *const pClusterW) const
Whether a trio of clusters are consistent with representing projections of the same 3d trajectory.
void SelectInputClusters(const pandora::ClusterList &inputClusterList, pandora::ClusterList &selectedClusterList) const
Select a subset of input clusters for processing in this algorithm.
void FindOverlaps(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2, SimpleOverlapTensor &overlapTensor) const
Find cluster overlaps and record these in the overlap tensor.
bool IsOverlap(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2) const
Whether two clusters overlap convincingly in x.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
Read the algorithm settings.
float m_minXOverlapFraction
The min x overlap fraction required in order to id overlapping clusters.
float m_minVertexLongitudinalDistance
Vertex association check: min longitudinal distance cut.
void CalculateEffectiveSpan(const pandora::Cluster *const pCluster, const float xMin, const float xMax, float &xMinEff, float &xMaxEff) const
Calculate effective span for a given clsuter taking gaps into account.
CartesianVector class.
Cluster class.
Definition Cluster.h:31
bool IsAvailable() const
Whether the cluster is available to be added to a particle flow object.
Definition Cluster.h:582
unsigned int GetNCaloHits() const
Get the number of calo hits in the cluster.
Definition Cluster.h:484
void GetClusterSpanX(float &xmin, float &xmax) const
Get minimum and maximum X positions of the calo hits in this cluster.
Definition Cluster.cc:169
const DetectorGapList & GetDetectorGapList() const
Get the list of gaps in the active detector volume.
bool ShouldDisplayAlgorithmInfo() const
Whether to display algorithm information during processing.
ParticleFlowObject class.
static float GetParticleMass(const int pdgCode)
Get the mass of a particle type.
Definition PdgTable.h:205
static int GetParticleCharge(const int pdgCode)
Get the charge of a particle type.
Definition PdgTable.h:227
const Pandora & GetPandora() const
Get the associated pandora instance.
Definition Process.h:116
StatusCodeException class.
static StatusCode ReadVectorOfValues(const TiXmlHandle &xmlHandle, const std::string &xmlElementName, std::vector< T > &vector)
Read a vector of values from a (space separated) list in an xml element.
Definition XmlHelper.h:229
static StatusCode ReadValue(const TiXmlHandle &xmlHandle, const std::string &xmlElementName, T &t)
Read a value from an xml element.
Definition XmlHelper.h:136
HitType
Calorimeter hit type enum.
std::vector< const Cluster * > ClusterVector
MANAGED_CONTAINER< const Cluster * > ClusterList
std::vector< CartesianVector > CartesianPointVector
StatusCode
The StatusCode enum.
MANAGED_CONTAINER< const ParticleFlowObject * > PfoList