LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
DeltaRayIdentificationAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
13 
15 
16 using namespace pandora;
17 
18 namespace lar_content
19 {
20 
21 DeltaRayIdentificationAlgorithm::DeltaRayIdentificationAlgorithm() :
22  m_distanceForMatching(3.f),
23  m_minParentLengthSquared(10.f * 10.f),
24  m_maxDaughterLengthSquared(175.f * 175.f)
25 {
26 }
27 
28 //------------------------------------------------------------------------------------------------------------------------------------------
29 
31 {
32  PfoVector parentPfos, daughterPfos;
33  this->GetPfos(m_parentPfoListName, parentPfos);
34  this->GetPfos(m_daughterPfoListName, daughterPfos);
35 
36  if (parentPfos.empty())
37  {
38  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
39  std::cout << "DeltaRayIdentificationAlgorithm: pfo list " << m_parentPfoListName << " unavailable." << std::endl;
40  return STATUS_CODE_SUCCESS;
41  }
42 
43  // Build parent/daughter associations (currently using length and proximity)
44  PfoAssociationMap pfoAssociationMap;
45  this->BuildAssociationMap(parentPfos, daughterPfos, pfoAssociationMap);
46 
47  // Create the parent/daughter links
48  PfoList newDaughterPfoList;
49  this->BuildParentDaughterLinks(pfoAssociationMap, newDaughterPfoList);
50 
51  if (!newDaughterPfoList.empty())
52  {
53  PANDORA_THROW_RESULT_IF(
54  STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*this, m_parentPfoListName, m_daughterPfoListName, newDaughterPfoList));
55  }
56 
57  return STATUS_CODE_SUCCESS;
58 }
59 
60 //------------------------------------------------------------------------------------------------------------------------------------------
61 
62 void DeltaRayIdentificationAlgorithm::GetPfos(const std::string &inputPfoListName, PfoVector &outputPfoVector) const
63 {
64  const PfoList *pPfoList = NULL;
65  PANDORA_THROW_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, inputPfoListName, pPfoList));
66 
67  if (NULL == pPfoList)
68  return;
69 
70  outputPfoVector.insert(outputPfoVector.end(), pPfoList->begin(), pPfoList->end());
71  std::sort(outputPfoVector.begin(), outputPfoVector.end(), LArPfoHelper::SortByNHits);
72 }
73 
74 //------------------------------------------------------------------------------------------------------------------------------------------
75 
76 void DeltaRayIdentificationAlgorithm::BuildAssociationMap(const PfoVector &parentPfos, const PfoVector &daughterPfos, PfoAssociationMap &pfoAssociationMap) const
77 {
78  PfoSet parentPfoList, daughterPfoList;
79  parentPfoList.insert(parentPfos.begin(), parentPfos.end());
80  daughterPfoList.insert(daughterPfos.begin(), daughterPfos.end());
81 
82  PfoVector allPfos(parentPfos.begin(), parentPfos.end());
83  allPfos.insert(allPfos.end(), daughterPfos.begin(), daughterPfos.end());
84 
85  // Loop over possible daughter Pfos in primary list
86  for (PfoVector::const_iterator iter1 = parentPfos.begin(), iterEnd1 = parentPfos.end(); iter1 != iterEnd1; ++iter1)
87  {
88  const ParticleFlowObject *const pDaughterPfo = *iter1;
89 
90  // Find the best parent Pfo using combined list
91  const ParticleFlowObject *pBestParentPfo = NULL;
92  float bestDisplacement(std::numeric_limits<float>::max());
93 
94  for (PfoVector::const_iterator iter2 = allPfos.begin(), iterEnd2 = allPfos.end(); iter2 != iterEnd2; ++iter2)
95  {
96  const ParticleFlowObject *const pThisParentPfo = *iter2;
97  float thisDisplacement(std::numeric_limits<float>::max());
98 
99  if (pDaughterPfo == pThisParentPfo)
100  continue;
101 
102  if (!this->IsAssociated(pDaughterPfo, pThisParentPfo, thisDisplacement))
103  continue;
104 
105  if (thisDisplacement < bestDisplacement)
106  {
107  bestDisplacement = thisDisplacement;
108  pBestParentPfo = pThisParentPfo;
109  }
110  }
111 
112  if (!pBestParentPfo)
113  continue;
114 
115  // Case 1: candidate parent comes from primary list
116  if (pBestParentPfo->GetParentPfoList().empty())
117  {
118  // Check: parent shouldn't live in the secondary list
119  if (daughterPfoList.count(pBestParentPfo))
120  throw StatusCodeException(STATUS_CODE_FAILURE);
121 
122  pfoAssociationMap.insert(PfoAssociationMap::value_type(pDaughterPfo, pBestParentPfo));
123  }
124 
125  // Case 2: candidate parent comes from secondary list
126  else
127  {
128  // Check: parent shouldn't live in the primary list
129  if (parentPfoList.count(pBestParentPfo))
130  throw StatusCodeException(STATUS_CODE_FAILURE);
131 
132  // Check: there should only be one parent
133  if (pBestParentPfo->GetParentPfoList().size() != 1)
134  throw StatusCodeException(STATUS_CODE_FAILURE);
135 
136  // Check: get the new parent (and check there is no grand-parent)
137  PfoList::const_iterator pIter = pBestParentPfo->GetParentPfoList().begin();
138  const ParticleFlowObject *const pReplacementParentPfo = *pIter;
139  if (pReplacementParentPfo->GetParentPfoList().size() != 0)
140  throw StatusCodeException(STATUS_CODE_FAILURE);
141 
142  pfoAssociationMap.insert(PfoAssociationMap::value_type(pDaughterPfo, pReplacementParentPfo));
143  }
144  }
145 }
146 
147 //------------------------------------------------------------------------------------------------------------------------------------------
148 
150  const ParticleFlowObject *const pDaughterPfo, const ParticleFlowObject *const pParentPfo, float &displacement) const
151 {
152  displacement = std::numeric_limits<float>::max();
153 
154  if (pDaughterPfo == pParentPfo)
155  return false;
156 
157  const float daughterLengthSquared(LArPfoHelper::GetTwoDLengthSquared(pDaughterPfo));
158  const float parentLengthSquared(LArPfoHelper::GetTwoDLengthSquared(pParentPfo));
159 
160  if (daughterLengthSquared > m_maxDaughterLengthSquared || parentLengthSquared < m_minParentLengthSquared || daughterLengthSquared > 0.5 * parentLengthSquared)
161  return false;
162 
163  const float transitionLengthSquared(125.f);
164  const float displacementCut((daughterLengthSquared > transitionLengthSquared)
166  : m_distanceForMatching * (2.f - daughterLengthSquared / transitionLengthSquared));
167 
168  try
169  {
170  displacement = this->GetTwoDSeparation(pDaughterPfo, pParentPfo);
171  }
172  catch (StatusCodeException &statusCodeException)
173  {
174  if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
175  throw statusCodeException;
176  }
177 
178  if (displacement > displacementCut)
179  return false;
180 
181  return true;
182 }
183 
184 //------------------------------------------------------------------------------------------------------------------------------------------
185 
186 float DeltaRayIdentificationAlgorithm::GetTwoDSeparation(const ParticleFlowObject *const pDaughterPfo, const ParticleFlowObject *const pParentPfo) const
187 {
188  CartesianPointVector vertexVectorU, vertexVectorV, vertexVectorW;
189  this->GetTwoDVertices(pDaughterPfo, TPC_VIEW_U, vertexVectorU);
190  this->GetTwoDVertices(pDaughterPfo, TPC_VIEW_V, vertexVectorV);
191  this->GetTwoDVertices(pDaughterPfo, TPC_VIEW_W, vertexVectorW);
192 
193  ClusterList clusterListU, clusterListV, clusterListW;
194  LArPfoHelper::GetClusters(pParentPfo, TPC_VIEW_U, clusterListU);
195  LArPfoHelper::GetClusters(pParentPfo, TPC_VIEW_V, clusterListV);
196  LArPfoHelper::GetClusters(pParentPfo, TPC_VIEW_W, clusterListW);
197 
198  float sumViews(0.f);
199  float sumDisplacementSquared(0.f);
200 
201  if (!vertexVectorU.empty())
202  {
203  const float thisDisplacement(this->GetClosestDistance(vertexVectorU, clusterListU));
204  sumDisplacementSquared += thisDisplacement * thisDisplacement;
205  sumViews += 1.f;
206  }
207 
208  if (!vertexVectorV.empty())
209  {
210  const float thisDisplacement(this->GetClosestDistance(vertexVectorV, clusterListV));
211  sumDisplacementSquared += thisDisplacement * thisDisplacement;
212  sumViews += 1.f;
213  }
214 
215  if (!vertexVectorW.empty())
216  {
217  const float thisDisplacement(this->GetClosestDistance(vertexVectorW, clusterListW));
218  sumDisplacementSquared += thisDisplacement * thisDisplacement;
219  sumViews += 1.f;
220  }
221 
222  if (sumViews < std::numeric_limits<float>::epsilon())
223  throw StatusCodeException(STATUS_CODE_FAILURE);
224 
225  return std::sqrt(sumDisplacementSquared / sumViews);
226 }
227 
228 //------------------------------------------------------------------------------------------------------------------------------------------
229 
230 void DeltaRayIdentificationAlgorithm::GetTwoDVertices(const ParticleFlowObject *const pPfo, const HitType &hitType, CartesianPointVector &vertexVector) const
231 {
232  ClusterList clusterList;
233  LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
234 
235  for (ClusterList::const_iterator iter = clusterList.begin(), iterEnd = clusterList.end(); iter != iterEnd; ++iter)
236  {
237  const Cluster *const pCluster = *iter;
238 
239  CartesianVector firstCoordinate(0.f, 0.f, 0.f), secondCoordinate(0.f, 0.f, 0.f);
240  LArClusterHelper::GetExtremalCoordinates(pCluster, firstCoordinate, secondCoordinate);
241 
242  vertexVector.push_back(firstCoordinate);
243  vertexVector.push_back(secondCoordinate);
244  }
245 }
246 
247 //------------------------------------------------------------------------------------------------------------------------------------------
248 
249 float DeltaRayIdentificationAlgorithm::GetClosestDistance(const CartesianPointVector &vertexVector, const ClusterList &clusterList) const
250 {
251  if (vertexVector.empty() || clusterList.empty())
252  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
253 
254  float bestDisplacement(std::numeric_limits<float>::max());
255 
256  for (CartesianPointVector::const_iterator iter1 = vertexVector.begin(), iterEnd1 = vertexVector.end(); iter1 != iterEnd1; ++iter1)
257  {
258  const CartesianVector &thisVertex = *iter1;
259 
260  for (ClusterList::const_iterator iter2 = clusterList.begin(), iterEnd2 = clusterList.end(); iter2 != iterEnd2; ++iter2)
261  {
262  const Cluster *const pCluster = *iter2;
263  const float thisDisplacement(LArClusterHelper::GetClosestDistance(thisVertex, pCluster));
264 
265  if (thisDisplacement < bestDisplacement)
266  bestDisplacement = thisDisplacement;
267  }
268  }
269 
270  return bestDisplacement;
271 }
272 
273 //------------------------------------------------------------------------------------------------------------------------------------------
274 
275 void DeltaRayIdentificationAlgorithm::BuildParentDaughterLinks(const PfoAssociationMap &pfoAssociationMap, PfoList &daughterPfoList) const
276 {
277  PfoList pfoList;
278  for (const auto &mapEntry : pfoAssociationMap)
279  pfoList.push_back(mapEntry.first);
280  pfoList.sort(LArPfoHelper::SortByNHits);
281 
282  for (const ParticleFlowObject *const pDaughterPfo : pfoList)
283  {
284  const ParticleFlowObject *const pParentPfo(this->GetParent(pfoAssociationMap, pDaughterPfo));
285 
286  if (!pParentPfo)
287  throw StatusCodeException(STATUS_CODE_FAILURE);
288 
289  if (!LArPfoHelper::IsTrack(pParentPfo))
290  continue;
291 
292  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SetPfoParentDaughterRelationship(*this, pParentPfo, pDaughterPfo));
293 
294  PandoraContentApi::ParticleFlowObject::Metadata metadata;
295  metadata.m_particleId = E_MINUS;
296  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::AlterMetadata(*this, pDaughterPfo, metadata));
297 
298  daughterPfoList.push_back(pDaughterPfo);
299  }
300 }
301 
302 //------------------------------------------------------------------------------------------------------------------------------------------
303 
304 const ParticleFlowObject *DeltaRayIdentificationAlgorithm::GetParent(const PfoAssociationMap &pfoAssociationMap, const ParticleFlowObject *const pPfo) const
305 {
306  const ParticleFlowObject *pParentPfo = nullptr;
307  const ParticleFlowObject *pDaughterPfo = pPfo;
308 
309  while (1)
310  {
311  PfoAssociationMap::const_iterator iter = pfoAssociationMap.find(pDaughterPfo);
312  if (pfoAssociationMap.end() == iter)
313  break;
314 
315  pParentPfo = iter->second;
316  pDaughterPfo = pParentPfo;
317  }
318 
319  return pParentPfo;
320 }
321 
322 //------------------------------------------------------------------------------------------------------------------------------------------
323 
324 StatusCode DeltaRayIdentificationAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
325 {
326  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "ParentPfoListName", m_parentPfoListName));
327  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "DaughterPfoListName", m_daughterPfoListName));
328 
329  PANDORA_RETURN_RESULT_IF_AND_IF(
330  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DistanceForMatching", m_distanceForMatching));
331 
332  float minParentLength = std::sqrt(m_minParentLengthSquared);
333  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinParentLength", minParentLength));
334  m_minParentLengthSquared = minParentLength * minParentLength;
335 
336  float maxDaughterLength = std::sqrt(m_maxDaughterLengthSquared);
337  PANDORA_RETURN_RESULT_IF_AND_IF(
338  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxDaughterLength", maxDaughterLength));
339  m_maxDaughterLengthSquared = maxDaughterLength * maxDaughterLength;
340 
341  return STATUS_CODE_SUCCESS;
342 }
343 
344 } // namespace lar_content
static bool SortByNHits(const pandora::ParticleFlowObject *const pLhs, const pandora::ParticleFlowObject *const pRhs)
Sort pfos by number of constituent hits.
Header file for the pfo helper class.
float m_distanceForMatching
Maximum allowed distance of delta ray from parent cosmic ray.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
static void GetClusters(const pandora::PfoList &pfoList, const pandora::HitType &hitType, pandora::ClusterList &clusterList)
Get a list of clusters of a particular hit type from a list of pfos.
void GetPfos(const std::string &inputPfoListName, pandora::PfoVector &outputPfoVector) const
Get the vector of Pfos, given the input list name.
const pandora::ParticleFlowObject * GetParent(const PfoAssociationMap &pfoAssociationMap, const pandora::ParticleFlowObject *const pPfo) const
For a given daughter, follow the parent/daughter links to find the overall parent.
float GetTwoDSeparation(const pandora::ParticleFlowObject *const pDaughterPfo, const pandora::ParticleFlowObject *const pParentPfo) const
Calculate 2D separation between two Pfos.
intermediate_table::const_iterator const_iterator
static bool IsTrack(const pandora::ParticleFlowObject *const pPfo)
Return track flag based on Pfo Particle ID.
std::string m_parentPfoListName
The parent pfo list name.
static float GetTwoDLengthSquared(const pandora::ParticleFlowObject *const pPfo)
Calculate length of Pfo using 2D clusters.
float m_maxDaughterLengthSquared
Maximum allowed length of daughter delta ray.
TFile f
Definition: plotHisto.C:6
std::string m_daughterPfoListName
The daughter pfo list name.
Header file for the cluster helper class.
void GetTwoDVertices(const pandora::ParticleFlowObject *const pPfo, const pandora::HitType &hitType, pandora::CartesianPointVector &vertexVector) const
Calculate 2D separation between two Pfos.
bool IsAssociated(const pandora::ParticleFlowObject *const pDaughterPfo, const pandora::ParticleFlowObject *const pParentPfo, float &displacement) const
Determine if a given pair of Pfos have a parent/daughter association.
static void GetExtremalCoordinates(const pandora::ClusterList &clusterList, pandora::CartesianVector &innerCoordinate, pandora::CartesianVector &outerCoordinate)
Get positions of the two most distant calo hits in a list of cluster (ordered by Z) ...
HitType
Definition: HitType.h:12
void BuildParentDaughterLinks(const PfoAssociationMap &pfoAssociationMap, pandora::PfoList &outputPfoList) const
Build the parent/daughter links from the map of parent/daughter associations.
void BuildAssociationMap(const pandora::PfoVector &inputPfos, const pandora::PfoVector &outputPfos, PfoAssociationMap &pfoAssociationMap) const
Build parent/daughter associations between PFOs.
float m_minParentLengthSquared
Minimum allowed length of parent cosmic ray.
float GetClosestDistance(const pandora::CartesianPointVector &vertexVector, const pandora::ClusterList &clusterList) const
Calculate closest 2D separation between a set of vertices and a set of clusters.
Header file for the delta ray identification algorithm class.
std::unordered_map< const pandora::ParticleFlowObject *, const pandora::ParticleFlowObject * > PfoAssociationMap
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.