LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CosmicRayBaseMatchingAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
15 
16 using namespace pandora;
17 
18 namespace lar_content
19 {
20 
21 StatusCode CosmicRayBaseMatchingAlgorithm::Run()
22 {
23  // Get the available clusters for each view
24  ClusterVector availableClustersU, availableClustersV, availableClustersW;
25  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNameU, availableClustersU));
26  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNameV, availableClustersV));
27  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNameW, availableClustersW));
28 
29  // Select clean clusters in each view
30  ClusterVector cleanClustersU, cleanClustersV, cleanClustersW;
31  this->SelectCleanClusters(availableClustersU, cleanClustersU);
32  this->SelectCleanClusters(availableClustersV, cleanClustersV);
33  this->SelectCleanClusters(availableClustersW, cleanClustersW);
34 
35  // Build associations between pairs of views
36  ClusterAssociationMap matchedClusterUV, matchedClusterVW, matchedClusterWU;
37  this->MatchClusters(cleanClustersU, cleanClustersV, matchedClusterUV);
38  this->MatchClusters(cleanClustersV, cleanClustersW, matchedClusterVW);
39  this->MatchClusters(cleanClustersW, cleanClustersU, matchedClusterWU);
40 
41  // Build particles from associations
42  ParticleList particleList;
43  this->MatchThreeViews(matchedClusterUV, matchedClusterVW, matchedClusterWU, particleList);
44  this->MatchTwoViews(matchedClusterUV, matchedClusterVW, matchedClusterWU, particleList);
45  this->BuildParticles(particleList);
46 
47  return STATUS_CODE_SUCCESS;
48 }
49 
50 //------------------------------------------------------------------------------------------------------------------------------------------
51 
52 StatusCode CosmicRayBaseMatchingAlgorithm::GetAvailableClusters(const std::string inputClusterListName, ClusterVector &clusterVector) const
53 {
54  const ClusterList *pClusterList = NULL;
55  PANDORA_RETURN_RESULT_IF_AND_IF(
56  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, inputClusterListName, pClusterList))
57 
58  if (!pClusterList || pClusterList->empty())
59  {
60  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
61  std::cout << "CosmicRayBaseMatchingAlgorithm: unable to find cluster list " << inputClusterListName << std::endl;
62 
63  return STATUS_CODE_SUCCESS;
64  }
65 
66  for (const Cluster *const pCluster : *pClusterList)
67  {
68  if (!pCluster->IsAvailable())
69  continue;
70 
71  clusterVector.push_back(pCluster);
72  }
73 
74  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
75 
76  return STATUS_CODE_SUCCESS;
77 }
78 
79 //------------------------------------------------------------------------------------------------------------------------------------------
80 
81 void CosmicRayBaseMatchingAlgorithm::MatchClusters(
82  const ClusterVector &clusterVector1, const ClusterVector &clusterVector2, ClusterAssociationMap &matchedClusters12) const
83 {
84  // Check that there are input clusters from both views
85  if (clusterVector1.empty() || clusterVector2.empty())
86  return;
87 
88  const HitType hitType1(LArClusterHelper::GetClusterHitType(*clusterVector1.begin()));
89  const HitType hitType2(LArClusterHelper::GetClusterHitType(*clusterVector2.begin()));
90 
91  if (hitType1 == hitType2)
92  throw StatusCodeException(STATUS_CODE_FAILURE);
93 
94  for (const Cluster *const pCluster1 : clusterVector1)
95  {
96  for (const Cluster *const pCluster2 : clusterVector2)
97  {
98  if (this->MatchClusters(pCluster1, pCluster2))
99  {
100  UIntSet daughterVolumeIntersection;
101  LArGeometryHelper::GetCommonDaughterVolumes(pCluster1, pCluster2, daughterVolumeIntersection);
102 
103  if (!daughterVolumeIntersection.empty())
104  matchedClusters12[pCluster1].push_back(pCluster2);
105  }
106  }
107  }
108 }
109 
110 //------------------------------------------------------------------------------------------------------------------------------------------
111 
112 void CosmicRayBaseMatchingAlgorithm::MatchThreeViews(const ClusterAssociationMap &matchedClusters12,
113  const ClusterAssociationMap &matchedClusters23, const ClusterAssociationMap &matchedClusters31, ParticleList &matchedParticles) const
114 {
115  if (matchedClusters12.empty() || matchedClusters23.empty() || matchedClusters31.empty())
116  return;
117 
118  ParticleList candidateParticles;
119 
120  ClusterList clusterList1;
121  for (const auto &mapEntry : matchedClusters12)
122  clusterList1.push_back(mapEntry.first);
123  clusterList1.sort(LArClusterHelper::SortByNHits);
124 
125  for (const Cluster *const pCluster1 : clusterList1)
126  {
127  const ClusterList &clusterList2(matchedClusters12.at(pCluster1));
128 
129  for (const Cluster *const pCluster2 : clusterList2)
130  {
131  ClusterAssociationMap::const_iterator iter23 = matchedClusters23.find(pCluster2);
132 
133  if (matchedClusters23.end() == iter23)
134  continue;
135 
136  const ClusterList &clusterList3 = iter23->second;
137 
138  for (const Cluster *const pCluster3 : clusterList3)
139  {
140  ClusterAssociationMap::const_iterator iter31 = matchedClusters31.find(pCluster3);
141 
142  if (matchedClusters31.end() == iter31)
143  continue;
144 
145  if (iter31->second.end() == std::find(iter31->second.begin(), iter31->second.end(), pCluster1))
146  continue;
147 
148  const HitType hitType1(LArClusterHelper::GetClusterHitType(pCluster1));
149  const HitType hitType2(LArClusterHelper::GetClusterHitType(pCluster2));
150  const HitType hitType3(LArClusterHelper::GetClusterHitType(pCluster3));
151 
152  if (!this->CheckMatchedClusters3D(pCluster1, pCluster2, pCluster3))
153  continue;
154 
155  const Cluster *const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1
156  : (TPC_VIEW_U == hitType2) ? pCluster2
157  : (TPC_VIEW_U == hitType3) ? pCluster3
158  : NULL);
159  const Cluster *const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1
160  : (TPC_VIEW_V == hitType2) ? pCluster2
161  : (TPC_VIEW_V == hitType3) ? pCluster3
162  : NULL);
163  const Cluster *const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1
164  : (TPC_VIEW_W == hitType2) ? pCluster2
165  : (TPC_VIEW_W == hitType3) ? pCluster3
166  : NULL);
167 
168  candidateParticles.push_back(Particle(pClusterU, pClusterV, pClusterW));
169  }
170  }
171  }
172 
173  return this->ResolveAmbiguities(candidateParticles, matchedParticles);
174 }
175 
176 //------------------------------------------------------------------------------------------------------------------------------------------
177 
178 void CosmicRayBaseMatchingAlgorithm::MatchTwoViews(const ClusterAssociationMap &matchedClusters12,
179  const ClusterAssociationMap &matchedClusters23, const ClusterAssociationMap &matchedClusters31, ParticleList &matchedParticles) const
180 {
181  ParticleList candidateParticles;
182  this->MatchTwoViews(matchedClusters12, candidateParticles);
183  this->MatchTwoViews(matchedClusters23, candidateParticles);
184  this->MatchTwoViews(matchedClusters31, candidateParticles);
185 
186  return this->ResolveAmbiguities(candidateParticles, matchedParticles);
187 }
188 
189 //------------------------------------------------------------------------------------------------------------------------------------------
190 
191 void CosmicRayBaseMatchingAlgorithm::MatchTwoViews(const ClusterAssociationMap &matchedClusters12, ParticleList &matchedParticles) const
192 {
193  if (matchedClusters12.empty())
194  return;
195 
196  ClusterList clusterList1;
197  for (const auto &mapEntry : matchedClusters12)
198  clusterList1.push_back(mapEntry.first);
199  clusterList1.sort(LArClusterHelper::SortByNHits);
200 
201  for (const Cluster *const pCluster1 : clusterList1)
202  {
203  const ClusterList &clusterList2(matchedClusters12.at(pCluster1));
204 
205  for (const Cluster *const pCluster2 : clusterList2)
206  {
207  const HitType hitType1(LArClusterHelper::GetClusterHitType(pCluster1));
208  const HitType hitType2(LArClusterHelper::GetClusterHitType(pCluster2));
209 
210  const Cluster *const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 : NULL);
211  const Cluster *const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 : NULL);
212  const Cluster *const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 : NULL);
213 
214  matchedParticles.push_back(Particle(pClusterU, pClusterV, pClusterW));
215  }
216  }
217 }
218 
219 //------------------------------------------------------------------------------------------------------------------------------------------
220 
221 void CosmicRayBaseMatchingAlgorithm::ResolveAmbiguities(const ParticleList &candidateParticles, ParticleList &matchedParticles) const
222 {
223  for (const Particle &particle1 : candidateParticles)
224  {
225  bool isGoodMatch(true);
226 
227  for (const Particle &particle2 : candidateParticles)
228  {
229  const bool commonU(particle1.m_pClusterU == particle2.m_pClusterU);
230  const bool commonV(particle1.m_pClusterV == particle2.m_pClusterV);
231  const bool commonW(particle1.m_pClusterW == particle2.m_pClusterW);
232 
233  const bool ambiguousU(commonU && NULL != particle1.m_pClusterU);
234  const bool ambiguousV(commonV && NULL != particle1.m_pClusterV);
235  const bool ambiguousW(commonW && NULL != particle1.m_pClusterW);
236 
237  if (commonU && commonV && commonW)
238  continue;
239 
240  if (ambiguousU || ambiguousV || ambiguousW)
241  {
242  isGoodMatch = false;
243  break;
244  }
245  }
246 
247  if (isGoodMatch)
248  matchedParticles.push_back(particle1);
249  }
250 }
251 
252 //------------------------------------------------------------------------------------------------------------------------------------------
253 
254 void CosmicRayBaseMatchingAlgorithm::BuildParticles(const ParticleList &particleList)
255 {
256  if (particleList.empty())
257  return;
258 
259  const PfoList *pPfoList = NULL;
260  std::string pfoListName;
261  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pPfoList, pfoListName));
262 
263  for (const Particle &particle : particleList)
264  {
265  const Cluster *const pClusterU = particle.m_pClusterU;
266  const Cluster *const pClusterV = particle.m_pClusterV;
267  const Cluster *const pClusterW = particle.m_pClusterW;
268 
269  const bool isAvailableU((NULL != pClusterU) ? pClusterU->IsAvailable() : true);
270  const bool isAvailableV((NULL != pClusterV) ? pClusterV->IsAvailable() : true);
271  const bool isAvailableW((NULL != pClusterW) ? pClusterW->IsAvailable() : true);
272 
273  if (!(isAvailableU && isAvailableV && isAvailableW))
274  throw StatusCodeException(STATUS_CODE_FAILURE);
275 
276  PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
277  this->SetPfoParameters(particle, pfoParameters);
278 
279  if (pfoParameters.m_clusterList.empty())
280  throw StatusCodeException(STATUS_CODE_FAILURE);
281 
282  const ParticleFlowObject *pPfo(NULL);
283  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pPfo));
284  }
285 
286  if (!pPfoList->empty())
287  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*this, m_outputPfoListName));
288 }
289 
290 //------------------------------------------------------------------------------------------------------------------------------------------
291 
292 CosmicRayBaseMatchingAlgorithm::Particle::Particle(const Cluster *const pClusterU, const Cluster *const pClusterV, const Cluster *const pClusterW) :
293  m_pClusterU(pClusterU),
294  m_pClusterV(pClusterV),
295  m_pClusterW(pClusterW)
296 {
297  if (NULL == m_pClusterU && NULL == m_pClusterV && NULL == m_pClusterW)
298  throw StatusCodeException(STATUS_CODE_FAILURE);
299 
300  const HitType hitTypeU(NULL == m_pClusterU ? TPC_VIEW_U : LArClusterHelper::GetClusterHitType(m_pClusterU));
301  const HitType hitTypeV(NULL == m_pClusterV ? TPC_VIEW_V : LArClusterHelper::GetClusterHitType(m_pClusterV));
302  const HitType hitTypeW(NULL == m_pClusterW ? TPC_VIEW_W : LArClusterHelper::GetClusterHitType(m_pClusterW));
303 
304  if (!(TPC_VIEW_U == hitTypeU && TPC_VIEW_V == hitTypeV && TPC_VIEW_W == hitTypeW))
305  throw StatusCodeException(STATUS_CODE_FAILURE);
306 }
307 
308 //------------------------------------------------------------------------------------------------------------------------------------------
309 
310 StatusCode CosmicRayBaseMatchingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
311 {
312  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputPfoListName", m_outputPfoListName));
313  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "InputClusterListNameU", m_inputClusterListNameU));
314  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "InputClusterListNameV", m_inputClusterListNameV));
315  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "InputClusterListNameW", m_inputClusterListNameW));
316 
317  return STATUS_CODE_SUCCESS;
318 }
319 
320 } // namespace lar_content
std::string m_inputClusterListNameV
The name of the view V cluster list.
Header file for the cosmic ray base matching algorithm class.
intermediate_table::const_iterator const_iterator
const pandora::Cluster * m_pClusterV
Address of cluster in V view.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
Header file for the geometry helper class.
Header file for the cluster helper class.
const pandora::Cluster * m_pClusterW
Address of cluster in W view.
std::string m_outputPfoListName
The name of the output PFO list.
const pandora::Cluster * m_pClusterU
Address of cluster in U view.
HitType
Definition: HitType.h:12
std::string m_inputClusterListNameU
The name of the view U cluster list.
std::string m_inputClusterListNameW
The name of the view W cluster list.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
std::vector< art::Ptr< recob::Cluster > > ClusterVector
std::unordered_map< const pandora::Cluster *, pandora::ClusterList > ClusterAssociationMap