LArSoft  v09_93_00
Liquid Argon Software toolkit - https://larsoft.org/
HierarchyValidationAlgorithm.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 HierarchyValidationAlgorithm::HierarchyValidationAlgorithm() :
22  m_event{-1},
23  m_detector{"dune_fd_hd"},
24  m_writeEventTree{false},
25  m_writeMCTree{false},
26  m_foldToPrimaries{false},
27  m_foldDynamic{false},
29  m_validateEvent{false},
30  m_validateMC{false},
31  m_minPurity{0.8f},
32  m_minCompleteness{0.65f},
33  m_minRecoHits{30},
37 {
38 }
39 
40 //------------------------------------------------------------------------------------------------------------------------------------------
41 
43 {
44  if (m_writeEventTree)
45  {
46  PANDORA_MONITORING_API(SaveTree(this->GetPandora(), m_eventTreeName.c_str(), m_eventFileName.c_str(), "UPDATE"));
47  }
48  if (m_writeMCTree)
49  {
50  PANDORA_MONITORING_API(SaveTree(this->GetPandora(), m_MCTreeName.c_str(), m_MCFileName.c_str(), "UPDATE"));
51  }
52 }
53 
54 //------------------------------------------------------------------------------------------------------------------------------------------
55 
57 {
58  ++m_event;
59  const CaloHitList *pCaloHitList(nullptr);
60  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*this, m_caloHitListName, pCaloHitList));
61  const MCParticleList *pMCParticleList(nullptr);
62  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pMCParticleList));
63  const PfoList *pPfoList(nullptr);
64  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*this, m_pfoListName, pPfoList));
65 
68  foldParameters.m_foldToTier = true;
69  else if (m_foldDynamic)
70  foldParameters.m_foldDynamic = true;
71  else if (m_foldToLeadingShowers)
72  foldParameters.m_foldToLeadingShowers = true;
75  LArHierarchyHelper::MCHierarchy mcHierarchy(recoCriteria);
76  LArHierarchyHelper::FillMCHierarchy(*pMCParticleList, *pCaloHitList, foldParameters, mcHierarchy);
78  LArHierarchyHelper::FillRecoHierarchy(*pPfoList, foldParameters, recoHierarchy);
80  LArHierarchyHelper::MatchInfo matchInfo(mcHierarchy, recoHierarchy, quality);
82  matchInfo.Print(mcHierarchy);
83 
84 #ifdef MONITORING
85  if (m_validateEvent)
86  this->EventValidation(matchInfo);
87  if (m_validateMC)
88  this->MCValidation(matchInfo);
89 #endif
90 
91  return STATUS_CODE_SUCCESS;
92 }
93 
94 //------------------------------------------------------------------------------------------------------------------------------------------
95 
96 #ifdef MONITORING
97 void HierarchyValidationAlgorithm::EventValidation(const LArHierarchyHelper::MatchInfo &matchInfo) const
98 {
99  if (m_writeEventTree)
100  {
101  int interaction{0};
102  MCParticleList rootMCParticles;
103  matchInfo.GetRootMCParticles(rootMCParticles);
104 
105  const LArHierarchyHelper::RecoHierarchy &recoHierarchy{matchInfo.GetRecoHierarchy()};
106  PfoList rootPfos;
107  recoHierarchy.GetRootPfos(rootPfos);
108  std::map<const LArHierarchyHelper::RecoHierarchy::Node *, const ParticleFlowObject *> recoNodeToRootMap;
109  for (const ParticleFlowObject *const pRoot : rootPfos)
110  {
112  recoHierarchy.GetFlattenedNodes(pRoot, nodes);
113  for (const LArHierarchyHelper::RecoHierarchy::Node *pNode : nodes)
114  recoNodeToRootMap[pNode] = pRoot;
115  }
116 
117  std::set<const pandora::ParticleFlowObject *> matchedRecoSliceRoots;
118  for (const MCParticle *const pRoot : rootMCParticles)
119  {
120  const LArHierarchyHelper::MCMatchesVector &matches{matchInfo.GetMatches(pRoot)};
121 
122  MCParticleList primaries;
123  for (const LArHierarchyHelper::MCMatches &match : matches)
124  {
125  const LArHierarchyHelper::MCHierarchy::Node *pMCNode{match.GetMC()};
126  if (pMCNode->GetHierarchyTier() == 1)
127  {
128  const MCParticle *const pLeadingMC{pMCNode->GetLeadingMCParticle()};
129  primaries.emplace_back(pLeadingMC);
130  }
131  }
132  // Check primaries is not empty before proceeding
133  if (primaries.size() == 0)
134  continue;
135  primaries.sort(LArMCParticleHelper::SortByMomentum);
137 
138  const int isCC{descriptor.IsCC()};
139  const int isQE{descriptor.IsQE()};
140  const int isResonant{descriptor.IsResonant()};
141  const int isDIS{descriptor.IsDIS()};
142  const int isCoherent{descriptor.IsCoherent()};
143  const int isNuMu{descriptor.IsMuonNeutrino()};
144  const int isNuE{descriptor.IsElectronNeutrino()};
145  const int nPiZero{static_cast<int>(descriptor.GetNumPiZero())};
146  const int nPiPlus{static_cast<int>(descriptor.GetNumPiPlus())};
147  const int nPiMinus{static_cast<int>(descriptor.GetNumPiMinus())};
148  const int nPhotons{static_cast<int>(descriptor.GetNumPhotons())};
149  const int nProtons{static_cast<int>(descriptor.GetNumProtons())};
150 
151  std::set<const LArHierarchyHelper::MCHierarchy::Node *> trackNodeSet, showerNodeSet;
152  int nGoodTrackMatches{0}, nGoodShowerMatches{0};
153  int nGoodMatches{0}, nPoorMatches{0}, nUnmatched{0};
154  int nGoodTier1Matches{0}, nTier1Nodes{0};
155  int nGoodTier1TrackMatches{0}, nTier1TrackNodes{0};
156  int nGoodTier1ShowerMatches{0}, nTier1ShowerNodes{0};
157  int hasLeadingMuon{0}, hasLeadingElectron{0}, isLeadingLeptonCorrect{0};
158  for (const LArHierarchyHelper::MCMatches &mcMatch : matches)
159  {
160  const LArHierarchyHelper::MCHierarchy::Node *pNode{mcMatch.GetMC()};
161 
162  for (const LArHierarchyHelper::RecoHierarchy::Node *pReco : mcMatch.GetRecoMatches())
163  matchedRecoSliceRoots.insert(recoNodeToRootMap[pReco]);
164 
165  const int nReco{static_cast<int>(mcMatch.GetRecoMatches().size())};
166  const bool isQuality{mcMatch.IsQuality(matchInfo.GetQualityCuts())};
167  if (nReco == 1 && isQuality)
168  ++nGoodMatches;
169  else if (nReco >= 1)
170  ++nPoorMatches;
171  else
172  ++nUnmatched;
173  if (pNode->GetHierarchyTier() == 1)
174  {
175  ++nTier1Nodes;
176  if (nReco == 1 && isQuality)
177  ++nGoodTier1Matches;
178  }
179 
180  const int pdg{std::abs(pNode->GetParticleId())};
181  if (pNode->IsLeadingLepton())
182  {
183  if (pdg == MU_MINUS)
184  hasLeadingMuon = 1;
185  else if (pdg == E_MINUS)
186  hasLeadingElectron = 1;
187  isLeadingLeptonCorrect = nReco == 1 ? 1 : 0;
188  }
189 
190  if (pdg == PHOTON || pdg == E_MINUS)
191  {
192  showerNodeSet.insert(pNode);
193  if (nReco == 1 && isQuality)
194  {
195  ++nGoodShowerMatches;
196  if (pNode->GetHierarchyTier() == 1)
197  ++nGoodTier1ShowerMatches;
198  }
199  if (pNode->GetHierarchyTier() == 1)
200  ++nTier1ShowerNodes;
201  }
202  else
203  {
204  trackNodeSet.insert(pNode);
205  if (nReco == 1 && isQuality)
206  {
207  ++nGoodTrackMatches;
208  if (pNode->GetHierarchyTier() == 1)
209  ++nGoodTier1TrackMatches;
210  }
211  if (pNode->GetHierarchyTier() == 1)
212  ++nTier1TrackNodes;
213  }
214  }
215 
216  const int nNodes{static_cast<int>(matchInfo.GetNMCNodes(pRoot))};
217  const int nTrackNodes{static_cast<int>(trackNodeSet.size())}, nShowerNodes{static_cast<int>(showerNodeSet.size())};
218  const int nRecoSlices{static_cast<int>(matchedRecoSliceRoots.size())};
219  const CartesianVector &trueVertex{pRoot->GetVertex()};
220  const float trueVtxX{trueVertex.GetX()};
221  const float trueVtxY{trueVertex.GetY()};
222  const float trueVtxZ{trueVertex.GetZ()};
223  float recoVtxX{std::numeric_limits<float>::max()};
224  float recoVtxY{std::numeric_limits<float>::max()};
225  float recoVtxZ{std::numeric_limits<float>::max()};
226  float vtxDx{std::numeric_limits<float>::max()};
227  float vtxDy{std::numeric_limits<float>::max()};
228  float vtxDz{std::numeric_limits<float>::max()};
229  float vtxDr{std::numeric_limits<float>::max()};
230  const int isFiducial{LArVertexHelper::IsInFiducialVolume(this->GetPandora(), trueVertex, m_detector)};
231 
232  for (const ParticleFlowObject *pRootPfo : matchedRecoSliceRoots)
233  {
234  try
235  {
236  const CartesianVector &recoVertex{LArPfoHelper::GetVertex(pRootPfo)->GetPosition()};
237  const float recoVertexX{recoVertex.GetX()};
238  const float recoVertexY{recoVertex.GetY()};
239  const float recoVertexZ{recoVertex.GetZ()};
240  const float dx{recoVertexX - trueVtxX};
241  const float dy{recoVertexY - trueVtxY};
242  const float dz{recoVertexZ - trueVtxZ};
243  const float dr{std::sqrt(dx * dx + dy * dy + dz * dz)};
244  if (dr < vtxDr)
245  {
246  recoVtxX = recoVertexX;
247  recoVtxY = recoVertexY;
248  recoVtxZ = recoVertexZ;
249  vtxDx = dx;
250  vtxDy = dy;
251  vtxDz = dz;
252  vtxDr = dr;
253  }
254  }
255  catch (StatusCodeException &)
256  {
257  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
258  std::cout << "HierarchyValidationAlgorithm::EventValidation could not retrieve recoVertex" << std::endl;
259  }
260  }
261 
262  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "event", m_event));
263  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "interaction", interaction));
264  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nRecoSlices", nRecoSlices));
265  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "isCC", isCC));
266  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "isQE", isQE));
267  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "isResonant", isResonant));
268  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "isDIS", isDIS));
269  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "isCoherent", isCoherent));
270  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "isNuMu", isNuMu));
271  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "isNuE", isNuE));
272  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nPiZero", nPiZero));
273  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nPiPlus", nPiPlus));
274  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nPiMinus", nPiMinus));
275  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nPhotons", nPhotons));
276  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nProtons", nProtons));
277  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "isFiducial", isFiducial));
278  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "recoVtxX", recoVtxX));
279  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "recoVtxY", recoVtxY));
280  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "recoVtxZ", recoVtxZ));
281  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "trueVtxX", recoVtxX));
282  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "trueVtxY", recoVtxY));
283  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "trueVtxZ", recoVtxZ));
284  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "vtxDx", vtxDx));
285  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "vtxDy", vtxDy));
286  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "vtxDz", vtxDz));
287  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "vtxDr", vtxDr));
288  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nGoodMatches", nGoodMatches));
289  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nPoorMatches", nPoorMatches));
290  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nUnmatched", nUnmatched));
291  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nNodes", nNodes));
292  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nGoodTier1Matches", nGoodTier1Matches));
293  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nTier1Nodes", nTier1Nodes));
294  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nGoodTrackMatches", nGoodTrackMatches));
295  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nGoodShowerMatches", nGoodShowerMatches));
296  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nTrackNodes", nTrackNodes));
297  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nShowerNodes", nShowerNodes));
298  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nGoodTier1TrackMatches", nGoodTier1TrackMatches));
299  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nTier1TrackNodes", nTier1TrackNodes));
300  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nGoodTier1ShowerMatches", nGoodTier1ShowerMatches));
301  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "nTier1ShowerNodes", nTier1ShowerNodes));
302  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "hasLeadingMuon", hasLeadingMuon));
303  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "hasLeadingElectron", hasLeadingElectron));
304  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_eventTreeName.c_str(), "isLeadingLeptonCorrect", isLeadingLeptonCorrect));
305  PANDORA_MONITORING_API(FillTree(this->GetPandora(), m_eventTreeName.c_str()));
306  ++interaction;
307  }
308  }
309 }
310 
311 //------------------------------------------------------------------------------------------------------------------------------------------
312 
313 void HierarchyValidationAlgorithm::MCValidation(const LArHierarchyHelper::MatchInfo &matchInfo) const
314 {
315  if (m_writeMCTree)
316  {
317  int interaction{0};
318  MCParticleList rootMCParticles;
319  matchInfo.GetRootMCParticles(rootMCParticles);
320  PfoList rootPfos;
321  const LArHierarchyHelper::RecoHierarchy &recoHierarchy{matchInfo.GetRecoHierarchy()};
322  recoHierarchy.GetRootPfos(rootPfos);
323  std::map<const LArHierarchyHelper::RecoHierarchy::Node *, const ParticleFlowObject *> recoNodeToRootMap;
324  for (const ParticleFlowObject *const pRoot : rootPfos)
325  {
327  recoHierarchy.GetFlattenedNodes(pRoot, nodes);
328  for (const LArHierarchyHelper::RecoHierarchy::Node *pNode : nodes)
329  recoNodeToRootMap[pNode] = pRoot;
330  }
331 
332  for (const MCParticle *const pRoot : rootMCParticles)
333  {
334  MCParticleList primaries;
335  for (const LArHierarchyHelper::MCMatches &matches : matchInfo.GetMatches(pRoot))
336  {
337  const LArHierarchyHelper::MCHierarchy::Node *pMCNode{matches.GetMC()};
338  if (pMCNode->GetHierarchyTier() == 1)
339  {
340  const MCParticle *const pLeadingMC{pMCNode->GetLeadingMCParticle()};
341  primaries.emplace_back(pLeadingMC);
342  }
343  }
344  // Check primaries is not empty before proceeding
345  if (primaries.size() == 0)
346  continue;
347  primaries.sort(LArMCParticleHelper::SortByMomentum);
349 
350  for (const LArHierarchyHelper::MCMatches &matches : matchInfo.GetMatches(pRoot))
351  {
352  const LArHierarchyHelper::MCHierarchy::Node *pMCNode{matches.GetMC()};
353  const int isTestBeam{pMCNode->IsTestBeamParticle() ? 1 : 0};
354  const int isCosmicRay{!isTestBeam && pMCNode->IsCosmicRay() ? 1 : 0};
355  const int isNeutrinoInt{!(isTestBeam || isCosmicRay) ? 1 : 0};
356  const int mcId{pMCNode->GetId()};
357  const int pdg{pMCNode->GetParticleId()};
358  const int tier{pMCNode->GetHierarchyTier()};
359  const int mcHits{static_cast<int>(pMCNode->GetCaloHits().size())};
360  const int isLeadingLepton{pMCNode->IsLeadingLepton() ? 1 : 0};
361 
362  const MCParticle *const pLeadingMC{pMCNode->GetLeadingMCParticle()};
363  const MCParticleList &parentList{pLeadingMC->GetParentList()};
364  const int isElectron{std::abs(pLeadingMC->GetParticleId()) == E_MINUS ? 1 : 0};
365  const int hasMuonParent{parentList.size() == 1 && std::abs(parentList.front()->GetParticleId()) == MU_MINUS ? 1 : 0};
366  const int isMichel{isElectron && hasMuonParent && LArMCParticleHelper::IsDecay(pLeadingMC) ? 1 : 0};
367  const CartesianVector mcMomentum{pLeadingMC->GetMomentum()};
368  const float mcMtm{mcMomentum.GetMagnitude()};
369  const float mcMtmX{mcMomentum.GetX()};
370  const float mcMtmY{mcMomentum.GetY()};
371  const float mcMtmZ{mcMomentum.GetZ()};
372  const float mcEnergy{pLeadingMC->GetEnergy()};
373 
374  const LArHierarchyHelper::RecoHierarchy::NodeVector &nodeVector{matches.GetRecoMatches()};
375  const int nMatches{static_cast<int>(nodeVector.size())};
376  IntVector recoSliceIdVector, recoIdVector, nRecoHitsVector, nSharedHitsVector;
377  FloatVector purityVector, completenessVector;
378  FloatVector purityAdcVector, completenessAdcVector;
379  FloatVector purityVectorU, purityVectorV, purityVectorW, completenessVectorU, completenessVectorV, completenessVectorW;
380  FloatVector purityAdcVectorU, purityAdcVectorV, purityAdcVectorW, completenessAdcVectorU, completenessAdcVectorV, completenessAdcVectorW;
381  const CartesianVector &trueVertex{pLeadingMC->GetVertex()};
382  const float trueVtxX{trueVertex.GetX()};
383  const float trueVtxY{trueVertex.GetY()};
384  const float trueVtxZ{trueVertex.GetZ()};
385  float recoVtxX{std::numeric_limits<float>::max()};
386  float recoVtxY{std::numeric_limits<float>::max()};
387  float recoVtxZ{std::numeric_limits<float>::max()};
388  float vtxDx{std::numeric_limits<float>::max()};
389  float vtxDy{std::numeric_limits<float>::max()};
390  float vtxDz{std::numeric_limits<float>::max()};
391  float vtxDr{std::numeric_limits<float>::max()};
392 
393  const int isCC{descriptor.IsCC()};
394  const int isQE{descriptor.IsQE()};
395  const int isResonant{descriptor.IsResonant()};
396  const int isDIS{descriptor.IsDIS()};
397  const int isCoherent{descriptor.IsCoherent()};
398  const int isNuMu{descriptor.IsMuonNeutrino()};
399  const int isNuE{descriptor.IsElectronNeutrino()};
400  const int nPiZero{static_cast<int>(descriptor.GetNumPiZero())};
401  const int nPiPlus{static_cast<int>(descriptor.GetNumPiPlus())};
402  const int nPiMinus{static_cast<int>(descriptor.GetNumPiMinus())};
403  const int nPhotons{static_cast<int>(descriptor.GetNumPhotons())};
404  const int nProtons{static_cast<int>(descriptor.GetNumProtons())};
405 
406  for (const LArHierarchyHelper::RecoHierarchy::Node *pRecoNode : nodeVector)
407  {
408  const int sliceId{static_cast<int>(
409  std::distance(rootPfos.begin(), std::find(rootPfos.begin(), rootPfos.end(), recoNodeToRootMap[pRecoNode])))};
410  recoSliceIdVector.emplace_back(sliceId);
411  recoIdVector.emplace_back(pRecoNode->GetParticleId());
412  nRecoHitsVector.emplace_back(static_cast<int>(pRecoNode->GetCaloHits().size()));
413  nSharedHitsVector.emplace_back(static_cast<int>(matches.GetSharedHits(pRecoNode)));
414  purityVector.emplace_back(matches.GetPurity(pRecoNode));
415  completenessVector.emplace_back(matches.GetCompleteness(pRecoNode));
416  purityAdcVector.emplace_back(matches.GetPurity(pRecoNode, true));
417  completenessAdcVector.emplace_back(matches.GetCompleteness(pRecoNode, true));
418  purityVectorU.emplace_back(matches.GetPurity(pRecoNode, TPC_VIEW_U));
419  purityVectorV.emplace_back(matches.GetPurity(pRecoNode, TPC_VIEW_V));
420  purityVectorW.emplace_back(matches.GetPurity(pRecoNode, TPC_VIEW_W));
421  completenessVectorU.emplace_back(matches.GetCompleteness(pRecoNode, TPC_VIEW_U));
422  completenessVectorV.emplace_back(matches.GetCompleteness(pRecoNode, TPC_VIEW_V));
423  completenessVectorW.emplace_back(matches.GetCompleteness(pRecoNode, TPC_VIEW_W));
424  purityAdcVectorU.emplace_back(matches.GetPurity(pRecoNode, TPC_VIEW_U, true));
425  purityAdcVectorV.emplace_back(matches.GetPurity(pRecoNode, TPC_VIEW_V, true));
426  purityAdcVectorW.emplace_back(matches.GetPurity(pRecoNode, TPC_VIEW_W, true));
427  completenessAdcVectorU.emplace_back(matches.GetCompleteness(pRecoNode, TPC_VIEW_U, true));
428  completenessAdcVectorV.emplace_back(matches.GetCompleteness(pRecoNode, TPC_VIEW_V, true));
429  completenessAdcVectorW.emplace_back(matches.GetCompleteness(pRecoNode, TPC_VIEW_W, true));
430  if (nMatches > 0)
431  {
432  try
433  {
434  const CartesianVector recoVertex{LArPfoHelper::GetVertex(pRecoNode->GetLeadingPfo())->GetPosition()};
435  recoVtxX = recoVertex.GetX();
436  recoVtxY = recoVertex.GetY();
437  recoVtxZ = recoVertex.GetZ();
438  vtxDx = recoVtxX - trueVtxX;
439  vtxDy = recoVtxY - trueVtxY;
440  vtxDz = recoVtxZ - trueVtxZ;
441  vtxDr = std::sqrt(vtxDx * vtxDx + vtxDy * vtxDy + vtxDz * vtxDz);
442  }
443  catch (StatusCodeException &)
444  {
445  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
446  std::cout << "HierarchyValidationAlgorithm::MCValidation could not retrieve recoVertex" << std::endl;
447  }
448  }
449  }
450 
451  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "event", m_event));
452  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "interaction", interaction));
453  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "mcId", mcId));
454  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "mcPDG", pdg));
455  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "mcTier", tier));
456  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "mcNHits", mcHits));
457  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "mcMtm", mcMtm));
458  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "mcMtmX", mcMtmX));
459  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "mcMtmY", mcMtmY));
460  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "mcMtmZ", mcMtmZ));
461  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "mcEnergy", mcEnergy));
462  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isNuInteraction", isNeutrinoInt));
463  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isCosmicRay", isCosmicRay));
464  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isTestBeam", isTestBeam));
465  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isLeadingLepton", isLeadingLepton));
466  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isMichel", isMichel));
467  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "nMatches", nMatches));
468  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "recoSliceIdVector", &recoSliceIdVector));
469  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "recoIdVector", &recoIdVector));
470  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "nRecoHitsVector", &nRecoHitsVector));
471  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "nSharedHitsVector", &nSharedHitsVector));
472  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "purityVector", &purityVector));
473  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "completenessVector", &completenessVector));
474  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "purityAdcVector", &purityAdcVector));
475  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "completenessAdcVector", &completenessAdcVector));
476  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "purityVectorU", &purityVectorU));
477  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "purityVectorV", &purityVectorV));
478  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "purityVectorW", &purityVectorW));
479  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "completenessVectorU", &completenessVectorU));
480  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "completenessVectorV", &completenessVectorV));
481  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "completenessVectorW", &completenessVectorW));
482  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "purityAdcVectorU", &purityAdcVectorU));
483  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "purityAdcVectorV", &purityAdcVectorV));
484  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "purityAdcVectorW", &purityAdcVectorW));
485  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "completenessAdcVectorU", &completenessAdcVectorU));
486  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "completenessAdcVectorV", &completenessAdcVectorV));
487  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "completenessAdcVectorW", &completenessAdcVectorW));
488  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "recoVtxX", recoVtxX));
489  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "recoVtxY", recoVtxY));
490  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "recoVtxZ", recoVtxZ));
491  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "trueVtxX", trueVtxX));
492  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "trueVtxY", trueVtxY));
493  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "trueVtxZ", trueVtxZ));
494  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "vtxDx", vtxDx));
495  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "vtxDy", vtxDy));
496  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "vtxDz", vtxDz));
497  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "vtxDr", vtxDr));
498  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isCC", isCC));
499  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isQE", isQE));
500  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isResonant", isResonant));
501  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isDIS", isDIS));
502  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isCoherent", isCoherent));
503  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isNuMu", isNuMu));
504  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "isNuE", isNuE));
505  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "nPiZero", nPiZero));
506  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "nPiPlus", nPiPlus));
507  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "nPiMinus", nPiMinus));
508  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "nPhotons", nPhotons));
509  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_MCTreeName.c_str(), "nProtons", nProtons));
510  PANDORA_MONITORING_API(FillTree(this->GetPandora(), m_MCTreeName.c_str()));
511  }
512  ++interaction;
513  }
514  }
515 }
516 #endif
517 
518 //------------------------------------------------------------------------------------------------------------------------------------------
519 
520 StatusCode HierarchyValidationAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
521 {
522  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "CaloHitListName", m_caloHitListName));
523  if (m_caloHitListName.empty())
524  m_caloHitListName = "CaloHitList2D";
525  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "PfoListName", m_pfoListName));
526  if (m_pfoListName.empty())
527  m_pfoListName = "RecreatedPfos";
528 
529  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "Detector", m_detector));
530 
531  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ValidateEvent", m_validateEvent));
532  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ValidateMC", m_validateMC));
533 
534  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "WriteEventTree", m_writeEventTree));
535  if (m_writeEventTree)
536  {
537  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "EventFileName", m_eventFileName));
538  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "EventTreeName", m_eventTreeName));
539  }
540 
541  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "WriteMCTree", m_writeMCTree));
542  if (m_writeMCTree)
543  {
544  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "MCFileName", m_MCFileName));
545  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "MCTreeName", m_MCTreeName));
546  }
547 
548  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "FoldToPrimaries", m_foldToPrimaries));
549  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "FoldDynamic", m_foldDynamic));
550  PANDORA_RETURN_RESULT_IF_AND_IF(
551  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "FoldToLeadingShowers", m_foldToLeadingShowers));
552 
553  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinPurity", m_minPurity));
554  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCompleteness", m_minCompleteness));
555  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinRecoHits", m_minRecoHits));
556  PANDORA_RETURN_RESULT_IF_AND_IF(
557  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinRecoHitsPerView", m_minRecoHitsPerView));
558  PANDORA_RETURN_RESULT_IF_AND_IF(
559  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinRecoGoodViews", m_minRecoGoodViews));
560  PANDORA_RETURN_RESULT_IF_AND_IF(
561  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "RemoveRecoNeutrons", m_removeRecoNeutrons));
562 
563  return STATUS_CODE_SUCCESS;
564 }
565 
566 } // namespace lar_content
unsigned int m_minRecoGoodViews
Minimum number of reconstructed primary good views.
const MCMatchesVector & GetMatches(const pandora::MCParticle *const pRoot) const
Retrieve the vector of matches (this will include null matches - i.e. MC nodes with no corresponding ...
const RecoHierarchy & GetRecoHierarchy() const
Retrieve the reco hierarchy used for the matching.
bool m_foldToTier
Whether or not to apply folding based on particle tier.
bool m_foldToPrimaries
Whether or not to fold the hierarchy back to primary particles.
Header file for the interaction type helper class.
unsigned int m_minRecoHits
Minimum number of reconstructed primary good hits.
void GetRootMCParticles(pandora::MCParticleList &rootMCParticles) const
Retrieve the root MC particles of the interaction hierarchies.
bool m_foldToLeadingShowers
Whether or not to fold the hierarchy back to leading shower particles.
bool m_writeEventTree
Whether or not to output event validation information to a ROOT file.
static const pandora::Vertex * GetVertex(const pandora::ParticleFlowObject *const pPfo)
Get the pfo vertex.
bool m_writeMCTree
Whether or not to output MC validation information to a ROOT file.
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< int > IntVector
static void MatchHierarchies(MatchInfo &matchInfo)
Finds the matches between reconstructed and MC hierarchies.
float m_minPurity
Minimum purity to tag a node as being of good quality.
std::string m_MCFileName
The name of the MC ROOT file to write.
static bool SortByMomentum(const pandora::MCParticle *const pLhs, const pandora::MCParticle *const pRhs)
Sort mc particles by their momentum.
bool m_removeRecoNeutrons
Whether to remove reconstructed neutrons and their downstream particles.
static void FillRecoHierarchy(const pandora::PfoList &pfoList, const FoldingParameters &foldParameters, RecoHierarchy &hierarchy)
Fill a reconstructed hierarchy based on the specified folding criteria (see RecoHierarchy::FillHierar...
const pandora::MCParticle * GetLeadingMCParticle() const
Retrieve the leading MC particle associated with this node.
std::string m_eventTreeName
The name of the event ROOT tree to write.
float m_minCompleteness
Minimum completeness to tag a node as being of good quality.
static bool IsInFiducialVolume(const pandora::Pandora &pandora, const pandora::CartesianVector &vertex, const std::string &detector)
Determine if a vertex is within a detector&#39;s fiducial volume. This throws a STATUS_CODE_INVALID_PARAM...
std::string m_MCTreeName
The name of the MC ROOT tree to write.
static bool IsDecay(const pandora::MCParticle *const pMCParticle)
Check whether or not an MC particle comes from a decay process.
const QualityCuts & GetQualityCuts() const
Retrieve the quality cuts for matching.
std::string m_pfoListName
Name of input PFO list.
bool m_foldDynamic
Whether or not to use process and topological information to make folding decisions.
Header file for the vertex helper class.
bool m_foldDynamic
Whether or not to fold the hierarchy dynamically.
void Print(const MCHierarchy &mcHierarchy) const
Prints information about which reco nodes are matched to the MC nodes, information about hit sharing...
bool m_foldToLeadingShowers
Whether or not to fold shower children to the leading shower particle.
std::string m_caloHitListName
Name of input calo hit list.
std::vector< MCMatches > MCMatchesVector
std::string m_eventFileName
The name of the event ROOT file to write.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
void GetRootPfos(pandora::PfoList &rootPfos) const
Retrieve the root particle flow objects of the interaction hierarchies.
static void FillMCHierarchy(const pandora::MCParticleList &mcParticleList, const pandora::CaloHitList &caloHitList, const FoldingParameters &foldParameters, MCHierarchy &hierarchy)
Fill an MC hierarchy based on the specified folding criteria (see MCHierarchy::FillHierarchy for deta...
unsigned int GetNMCNodes(const pandora::MCParticle *const pRoot) const
Retrieve the number of MC nodes available to match.
unsigned int m_minRecoHitsPerView
Minimum number of reconstructed hits for a good view.
static InteractionDescriptor GetInteractionDescriptor(const pandora::MCParticleList &mcPrimaryList)
Get the interaction descriptor of an event.
bool m_validateMC
Whether to validate at the level of MC nodes.
bool IsCC() const
Whether or not the interaction is CC.
bool IsTestBeamParticle() const
Check if this is a test beam particle.
Header file for the hierarchy validation algorithm.
bool m_validateEvent
Whether to validate at the level of an event.