LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
TestBeamHierarchyEventValidationAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
14 
16 
17 #include <sstream>
18 
19 using namespace pandora;
20 
21 namespace lar_content
22 {
23 
24 TestBeamHierarchyEventValidationAlgorithm::TestBeamHierarchyEventValidationAlgorithm()
25 {
26 }
27 
28 //------------------------------------------------------------------------------------------------------------------------------------------
29 
30 TestBeamHierarchyEventValidationAlgorithm::~TestBeamHierarchyEventValidationAlgorithm()
31 {
32 }
33 
34 //------------------------------------------------------------------------------------------------------------------------------------------
35 
36 void TestBeamHierarchyEventValidationAlgorithm::FillValidationInfo(const MCParticleList *const pMCParticleList,
37  const CaloHitList *const pCaloHitList, const PfoList *const pPfoList, ValidationInfo &validationInfo) const
38 {
39  if (pMCParticleList && pCaloHitList)
40  {
41  LArMCParticleHelper::MCContributionMap targetMCParticleToHitsMap;
42  LArMCParticleHelper::SelectReconstructableTestBeamHierarchyMCParticles(
43  pMCParticleList, pCaloHitList, m_primaryParameters, LArMCParticleHelper::IsLeadingBeamParticle, targetMCParticleToHitsMap);
44  LArMCParticleHelper::SelectReconstructableMCParticles(
45  pMCParticleList, pCaloHitList, m_primaryParameters, LArMCParticleHelper::IsCosmicRay, targetMCParticleToHitsMap);
46 
47  LArMCParticleHelper::PrimaryParameters parameters(m_primaryParameters);
48  parameters.m_minPrimaryGoodHits = 0;
49  parameters.m_minHitsForGoodView = 0;
50  parameters.m_minHitSharingFraction = 0.f;
51  LArMCParticleHelper::MCContributionMap allMCParticleToHitsMap;
52  LArMCParticleHelper::SelectReconstructableTestBeamHierarchyMCParticles(
53  pMCParticleList, pCaloHitList, parameters, LArMCParticleHelper::IsLeadingBeamParticle, allMCParticleToHitsMap);
54  LArMCParticleHelper::SelectReconstructableMCParticles(
55  pMCParticleList, pCaloHitList, parameters, LArMCParticleHelper::IsCosmicRay, allMCParticleToHitsMap);
56 
57  validationInfo.SetTargetMCParticleToHitsMap(targetMCParticleToHitsMap);
58  validationInfo.SetAllMCParticleToHitsMap(allMCParticleToHitsMap);
59  }
60 
61  if (pPfoList)
62  {
63  PfoList allConnectedPfos;
64  LArPfoHelper::GetAllConnectedPfos(*pPfoList, allConnectedPfos);
65 
66  PfoList finalStatePfos;
67  for (const ParticleFlowObject *const pPfo : allConnectedPfos)
68  {
69  // ATTN: Is test beam only set for parent pfo, therefor add parent and daughters for that particle
70  if (LArPfoHelper::IsTestBeam(pPfo))
71  {
72  finalStatePfos.push_back(pPfo);
73  for (const ParticleFlowObject *const pDaughterPfo : pPfo->GetDaughterPfoList())
74  finalStatePfos.push_back(pDaughterPfo);
75  }
76  else if (pPfo->GetParentPfoList().empty())
77  {
78  finalStatePfos.push_back(pPfo);
79  }
80  }
81 
83  LArMCParticleHelper::GetTestBeamHierarchyPfoToReconstructable2DHitsMap(
84  finalStatePfos, validationInfo.GetAllMCParticleToHitsMap(), pfoToHitsMap, m_primaryParameters.m_foldBackHierarchy);
85  validationInfo.SetPfoToHitsMap(pfoToHitsMap);
86  }
87 
90  LArMCParticleHelper::GetPfoMCParticleHitSharingMaps(
91  validationInfo.GetPfoToHitsMap(), {validationInfo.GetAllMCParticleToHitsMap()}, pfoToMCHitSharingMap, mcToPfoHitSharingMap);
92  validationInfo.SetMCToPfoHitSharingMap(mcToPfoHitSharingMap);
93 
94  LArMCParticleHelper::MCParticleToPfoHitSharingMap interpretedMCToPfoHitSharingMap;
95  this->InterpretMatching(validationInfo, interpretedMCToPfoHitSharingMap);
96  validationInfo.SetInterpretedMCToPfoHitSharingMap(interpretedMCToPfoHitSharingMap);
97 }
98 
99 //------------------------------------------------------------------------------------------------------------------------------------------
100 
101 void TestBeamHierarchyEventValidationAlgorithm::ProcessOutput(
102  const ValidationInfo &validationInfo, const bool useInterpretedMatching, const bool printToScreen, const bool fillTree) const
103 {
104  static int eventNumber{-1};
105  ++eventNumber;
106  if (printToScreen && useInterpretedMatching)
107  std::cout << "---EVENT-" << eventNumber << "-INTERPRETED-MATCHING-OUTPUT----------------------------------------------------------"
108  << std::endl;
109  else if (printToScreen)
110  std::cout << "---EVENT-" << eventNumber << "-RAW-MATCHING-OUTPUT------------------------------------------------------------------"
111  << std::endl;
112 
113  PfoVector primaryPfoVector;
114  LArMonitoringHelper::GetOrderedPfoVector(validationInfo.GetPfoToHitsMap(), primaryPfoVector);
115 
116  // Test Beam Hierarchy Validation Pfo Bookkeeping
117  int pfoIndex(0), testBeamPfoIndex(0);
118  PfoToIdMap pfoToIdMap, testBeamPfoToIdMap;
119 
120  for (const Pfo *const pPrimaryPfo : primaryPfoVector)
121  {
122  pfoToIdMap.insert(PfoToIdMap::value_type(pPrimaryPfo, ++pfoIndex));
123  const Pfo *const pRecoTestBeam(LArPfoHelper::IsTestBeamFinalState(pPrimaryPfo) ? LArPfoHelper::GetParentPfo(pPrimaryPfo) : nullptr);
124 
125  if (pRecoTestBeam && !testBeamPfoToIdMap.count(pRecoTestBeam))
126  testBeamPfoToIdMap.insert(PfoToIdMap::value_type(pRecoTestBeam, ++testBeamPfoIndex));
127  }
128 
129  const LArMCParticleHelper::MCParticleToPfoHitSharingMap &mcToPfoHitSharingMap(
130  useInterpretedMatching ? validationInfo.GetInterpretedMCToPfoHitSharingMap() : validationInfo.GetMCToPfoHitSharingMap());
131 
132  // Test Beam Hierarchy Validation MCParticle Bookkeeping
133  MCParticleVector mcPrimaryVector;
134  LArMonitoringHelper::GetOrderedMCParticleVector({validationInfo.GetTargetMCParticleToHitsMap()}, mcPrimaryVector);
135  LArMCParticleHelper::MCParticleIntMap triggeredToLeading, triggeredToLeadingCounter;
136 
137  // ATTN: At this stage the mcPrimaryVector is ordered from neutrinos, beam and then cosmics. Here we extract the beam and reorder
138  // to ensure the order follows primary parent beam 1, daughter 1 of beam 1, daughter 2 of beam 1, ..., primary parent beam 2,
139  // daughter 1 of beam 2, etc... as expected by downstream logic
140  MCParticleVector mcPrimaryVectorCopy(mcPrimaryVector), triggeredBeamParticles;
141  LArMCParticleHelper::MCRelationMap leadingToTriggeredMap;
142  mcPrimaryVector.clear();
143 
144  for (const MCParticle *const pMCPrimary : mcPrimaryVectorCopy)
145  {
146  if (LArMCParticleHelper::IsLeadingBeamParticle(pMCPrimary))
147  {
148  const MCParticle *const pParentMCParticle(LArMCParticleHelper::GetParentMCParticle(pMCPrimary));
149  leadingToTriggeredMap.insert(LArMCParticleHelper::MCRelationMap::value_type(pMCPrimary, pParentMCParticle));
150 
151  if (std::find(triggeredBeamParticles.begin(), triggeredBeamParticles.end(), pParentMCParticle) == triggeredBeamParticles.end())
152  triggeredBeamParticles.push_back(pParentMCParticle);
153  }
154  else
155  {
156  mcPrimaryVector.push_back(pMCPrimary);
157  }
158  }
159 
160  for (const MCParticle *const pMCParent : triggeredBeamParticles)
161  {
162  // Parent appears first
163  mcPrimaryVector.push_back(pMCParent);
164  triggeredToLeading.insert(LArMCParticleHelper::MCParticleIntMap::value_type(pMCParent, 1));
165  triggeredToLeadingCounter.insert(LArMCParticleHelper::MCParticleIntMap::value_type(pMCParent, 0));
166 
167  for (const auto iter : leadingToTriggeredMap)
168  {
169  // Followed by daughters, veto parent <-> parent matche
170  if (iter.second == pMCParent && iter.first != pMCParent)
171  {
172  mcPrimaryVector.push_back(iter.first);
173  triggeredToLeading.at(pMCParent)++;
174  }
175  }
176  }
177 
178  PfoSet recoTestBeamHierarchies;
179  MCParticleList associatedMCPrimaries;
180 
181  int nCorrectTBHierarchy(0), nTotalTBHierarchy(0), nCorrectCR(0), nTotalCR(0);
182  int nFakeTBHierarchy(0), nFakeCR(0), nSplitTBHierarchy(0), nSplitCR(0), nLost(0), mcPrimaryIndex(0), nTargetMatches(0),
183  nTargetTBHierarchyMatches(0);
184  int nTargetCRMatches(0), nTargetGoodTBHierarchyMatches(0), nTargetTBHierarchySplits(0), nTargetTBHierarchyLosses(0);
185  IntVector mcPrimaryId, mcPrimaryPdg, mcPrimaryTier, nMCHitsTotal, nMCHitsU, nMCHitsV, nMCHitsW;
186  FloatVector mcPrimaryE, mcPrimaryPX, mcPrimaryPY, mcPrimaryPZ;
187  FloatVector mcPrimaryVtxX, mcPrimaryVtxY, mcPrimaryVtxZ, mcPrimaryEndX, mcPrimaryEndY, mcPrimaryEndZ;
188  IntVector nPrimaryMatchedPfos, nPrimaryMatchedTBHierarchyPfos, nPrimaryMatchedCRPfos;
189  IntVector bestMatchPfoId, bestMatchPfoPdg, bestMatchPfoTier, bestMatchPfoIsTestBeam, bestMatchPfoIsTestBeamHierarchy;
190  IntVector bestMatchPfoRecoTBId, bestMatchPfoNHitsTotal, bestMatchPfoNHitsU, bestMatchPfoNHitsV, bestMatchPfoNHitsW;
191  IntVector bestMatchPfoNSharedHitsTotal, bestMatchPfoNSharedHitsU, bestMatchPfoNSharedHitsV, bestMatchPfoNSharedHitsW;
192  FloatVector bestMatchPfoX0;
193 
194  std::stringstream targetSS;
195  const std::string name("TB");
196 
197  for (const MCParticle *const pMCPrimary : mcPrimaryVector)
198  {
199  const bool hasMatch(mcToPfoHitSharingMap.count(pMCPrimary) && !mcToPfoHitSharingMap.at(pMCPrimary).empty());
200  const bool isTargetPrimary(validationInfo.GetTargetMCParticleToHitsMap().count(pMCPrimary));
201 
202  if (!isTargetPrimary)
203  continue;
204 
205  // Parent in hierarchy needed even if no match
206  const bool hasVisibleTargets(
207  (!triggeredToLeading.empty() && LArMCParticleHelper::IsBeamParticle(LArMCParticleHelper::GetParentMCParticle(pMCPrimary)))
208  ? triggeredToLeading.at(LArMCParticleHelper::GetParentMCParticle(pMCPrimary)) != 1
209  : false);
210 
211  if (!hasMatch && !hasVisibleTargets)
212  continue;
213 
214  associatedMCPrimaries.push_back(pMCPrimary);
215  const int nTargetPrimaries(associatedMCPrimaries.size());
216  ++mcPrimaryIndex;
217  const CaloHitList &mcPrimaryHitList(validationInfo.GetAllMCParticleToHitsMap().at(pMCPrimary));
218 
219  const int mcNuanceCode(LArMCParticleHelper::GetNuanceCode(LArMCParticleHelper::GetParentMCParticle(pMCPrimary)));
220  const int isBeamParticle(LArMCParticleHelper::IsBeamParticle(pMCPrimary));
221 
222  // Leading beam particle is the primary beam particle or a daughter of that particle
223  const int isLeadingBeamParticle(LArMCParticleHelper::IsLeadingBeamParticle(pMCPrimary));
224  const int isCosmicRay(LArMCParticleHelper::IsCosmicRay(pMCPrimary));
225 
226  // Tier (0) : Primary, (1) : Daughter, (2) : Granddaughter etc... Note tier increases for both visible and invisible particles
227  const int mcHierarchyTier(LArMCParticleHelper::GetHierarchyTier(pMCPrimary));
228 
229  // Identify the number of matched leading particles and flag whether last particle in hierarchy is being considered
230  bool isLastTestBeamLeading(false);
231  if (isLeadingBeamParticle)
232  {
233  triggeredToLeadingCounter.at(LArMCParticleHelper::GetParentMCParticle(pMCPrimary))++;
234  const int nHierarchyLeading(triggeredToLeadingCounter.at(LArMCParticleHelper::GetParentMCParticle(pMCPrimary)));
235  isLastTestBeamLeading = (nHierarchyLeading == triggeredToLeading.at(LArMCParticleHelper::GetParentMCParticle(pMCPrimary)));
236  }
237 
238 #ifdef MONITORING
239  const CartesianVector &targetVertex(LArMCParticleHelper::GetParentMCParticle(pMCPrimary)->GetVertex());
240  const float targetVertexX(targetVertex.GetX()), targetVertexY(targetVertex.GetY()), targetVertexZ(targetVertex.GetZ());
241 #endif
242 
243  for (int tier = 0; tier < mcHierarchyTier; tier++)
244  targetSS << " -> ";
245 
246  targetSS << (!isTargetPrimary ? "(Non target) " : "") << "PrimaryId " << mcPrimaryIndex << ", TB " << isBeamParticle << ", TB Hierarchy "
247  << isLeadingBeamParticle << ", CR " << isCosmicRay << ", MCPDG " << pMCPrimary->GetParticleId() << ", Tier " << mcHierarchyTier
248  << ", Energy " << pMCPrimary->GetEnergy() << ", Dist. " << (pMCPrimary->GetEndpoint() - pMCPrimary->GetVertex()).GetMagnitude()
249  << ", nMCHits " << mcPrimaryHitList.size() << " (" << LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, mcPrimaryHitList)
250  << ", " << LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, mcPrimaryHitList) << ", "
251  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, mcPrimaryHitList) << ")" << std::endl;
252 
253  mcPrimaryId.push_back(mcPrimaryIndex);
254  mcPrimaryPdg.push_back(pMCPrimary->GetParticleId());
255  mcPrimaryTier.push_back(mcHierarchyTier);
256  mcPrimaryE.push_back(pMCPrimary->GetEnergy());
257  mcPrimaryPX.push_back(pMCPrimary->GetMomentum().GetX());
258  mcPrimaryPY.push_back(pMCPrimary->GetMomentum().GetY());
259  mcPrimaryPZ.push_back(pMCPrimary->GetMomentum().GetZ());
260  mcPrimaryVtxX.push_back(pMCPrimary->GetVertex().GetX());
261  mcPrimaryVtxY.push_back(pMCPrimary->GetVertex().GetY());
262  mcPrimaryVtxZ.push_back(pMCPrimary->GetVertex().GetZ());
263  mcPrimaryEndX.push_back(pMCPrimary->GetEndpoint().GetX());
264  mcPrimaryEndY.push_back(pMCPrimary->GetEndpoint().GetY());
265  mcPrimaryEndZ.push_back(pMCPrimary->GetEndpoint().GetZ());
266  nMCHitsTotal.push_back(mcPrimaryHitList.size());
267  nMCHitsU.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, mcPrimaryHitList));
268  nMCHitsV.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, mcPrimaryHitList));
269  nMCHitsW.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, mcPrimaryHitList));
270 
271  int matchIndex(0), nPrimaryMatches(0), nPrimaryTBHierarchyMatches(0), nPrimaryCRMatches(0), nPrimaryGoodTBHierarchyMatches(0),
272  nPrimaryTBHierarchySplits(0);
273 #ifdef MONITORING
274  float recoVertexX(std::numeric_limits<float>::max()), recoVertexY(std::numeric_limits<float>::max()),
275  recoVertexZ(std::numeric_limits<float>::max());
276 #endif
277  for (const LArMCParticleHelper::PfoCaloHitListPair &pfoToSharedHits : mcToPfoHitSharingMap.at(pMCPrimary))
278  {
279  const CaloHitList &sharedHitList(pfoToSharedHits.second);
280  const CaloHitList &pfoHitList(validationInfo.GetPfoToHitsMap().at(pfoToSharedHits.first));
281 
282  const bool isRecoTestBeam(LArPfoHelper::IsTestBeam(pfoToSharedHits.first));
283  const bool isRecoTestBeamHierarchy(LArPfoHelper::IsTestBeamFinalState(pfoToSharedHits.first));
284  const bool isGoodMatch(this->IsGoodMatch(mcPrimaryHitList, pfoHitList, sharedHitList));
285 
286  // Tier (0) : Primary, (1) : Daughter, (2) : Granddaughter etc... Note that the tier only increases for visible particle
287  const int pfoHierarchyTier(LArPfoHelper::GetHierarchyTier(pfoToSharedHits.first));
288  const int pfoId(pfoToIdMap.at(pfoToSharedHits.first));
289  const int recoTBId(
290  isRecoTestBeam || isRecoTestBeamHierarchy ? testBeamPfoToIdMap.at(LArPfoHelper::GetParentPfo(pfoToSharedHits.first)) : -1);
291 
292  if (0 == matchIndex++)
293  {
294  bestMatchPfoId.push_back(pfoId);
295  bestMatchPfoPdg.push_back(pfoToSharedHits.first->GetParticleId());
296  bestMatchPfoTier.push_back(pfoHierarchyTier);
297  bestMatchPfoIsTestBeam.push_back(isRecoTestBeam ? 1 : 0);
298  bestMatchPfoIsTestBeamHierarchy.push_back(isRecoTestBeamHierarchy ? 1 : 0);
299  bestMatchPfoRecoTBId.push_back(recoTBId);
300  bestMatchPfoNHitsTotal.push_back(pfoHitList.size());
301  bestMatchPfoNHitsU.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, pfoHitList));
302  bestMatchPfoNHitsV.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, pfoHitList));
303  bestMatchPfoNHitsW.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, pfoHitList));
304  bestMatchPfoNSharedHitsTotal.push_back(sharedHitList.size());
305  bestMatchPfoNSharedHitsU.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, sharedHitList));
306  bestMatchPfoNSharedHitsV.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, sharedHitList));
307  bestMatchPfoNSharedHitsW.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, sharedHitList));
308  bestMatchPfoX0.push_back(pfoToSharedHits.first->GetPropertiesMap().count("X0") ? pfoToSharedHits.first->GetPropertiesMap().at("X0")
309  : std::numeric_limits<float>::max());
310 #ifdef MONITORING
311  try
312  {
313  const Vertex *const pRecoVertex(isRecoTestBeamHierarchy
314  ? LArPfoHelper::GetTestBeamInteractionVertex(LArPfoHelper::GetParentPfo(pfoToSharedHits.first))
315  : LArPfoHelper::GetVertex(pfoToSharedHits.first));
316  recoVertexX = pRecoVertex->GetPosition().GetX();
317  recoVertexY = pRecoVertex->GetPosition().GetY();
318  recoVertexZ = pRecoVertex->GetPosition().GetZ();
319  }
320  catch (const StatusCodeException &)
321  {
322  }
323 #endif
324  }
325 
326  if (isGoodMatch)
327  ++nPrimaryMatches;
328 
329  // ATTN: In hierarchy mode let TBHierarchyMatches become effective TBHierarchyMatches and treat the same
330  if (isRecoTestBeamHierarchy && isGoodMatch)
331  ++nPrimaryTBHierarchyMatches;
332  if (!isRecoTestBeamHierarchy && isGoodMatch)
333  ++nPrimaryCRMatches;
334 
335  if (isRecoTestBeamHierarchy)
336  {
337  // Account for splitting of test beam particle into separate reconstructed primary pfos
338  const Pfo *const pRecoTB(LArPfoHelper::GetParentPfo(pfoToSharedHits.first));
339  const bool isSplitRecoTBHierarchy(!recoTestBeamHierarchies.empty() && !recoTestBeamHierarchies.count(pRecoTB));
340  if (!isSplitRecoTBHierarchy && isGoodMatch)
341  ++nPrimaryGoodTBHierarchyMatches;
342  if (isSplitRecoTBHierarchy && isLeadingBeamParticle && isGoodMatch)
343  ++nPrimaryTBHierarchySplits;
344  recoTestBeamHierarchies.insert(pRecoTB);
345  }
346 
347  for (int tier = 0; tier < mcHierarchyTier; tier++)
348  targetSS << " ";
349 
350  targetSS << "-" << (!isGoodMatch ? "(Below threshold) " : "") << "MatchedPfoId " << pfoId << ", TB " << isRecoTestBeam
351  << ", TB Hierarchy " << isRecoTestBeamHierarchy;
352  if (isRecoTestBeamHierarchy)
353  targetSS << " [TBId: " << recoTBId << "]";
354  targetSS << ", CR " << (!isRecoTestBeam && !isRecoTestBeamHierarchy) << ", PDG " << pfoToSharedHits.first->GetParticleId()
355  << ", Tier " << pfoHierarchyTier << ", nMatchedHits " << sharedHitList.size() << " ("
356  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, sharedHitList) << ", "
357  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, sharedHitList) << ", "
358  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, sharedHitList) << ")"
359  << ", nPfoHits " << pfoHitList.size() << " (" << LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, pfoHitList) << ", "
360  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, pfoHitList) << ", "
361  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, pfoHitList) << ")" << std::endl;
362  }
363 
364  if (mcToPfoHitSharingMap.at(pMCPrimary).empty())
365  {
366  for (int tier = 0; tier < mcHierarchyTier; tier++)
367  targetSS << " ";
368  targetSS << "-No matched Pfo" << std::endl;
369  bestMatchPfoId.push_back(-1);
370  bestMatchPfoPdg.push_back(0);
371  bestMatchPfoTier.push_back(-1);
372  bestMatchPfoIsTestBeam.push_back(0);
373  bestMatchPfoIsTestBeamHierarchy.push_back(0);
374  bestMatchPfoRecoTBId.push_back(-1);
375  bestMatchPfoNHitsTotal.push_back(0);
376  bestMatchPfoNHitsU.push_back(0);
377  bestMatchPfoNHitsV.push_back(0);
378  bestMatchPfoNHitsW.push_back(0);
379  bestMatchPfoNSharedHitsTotal.push_back(0);
380  bestMatchPfoNSharedHitsU.push_back(0);
381  bestMatchPfoNSharedHitsV.push_back(0);
382  bestMatchPfoNSharedHitsW.push_back(0);
383  bestMatchPfoX0.push_back(std::numeric_limits<float>::max());
384  }
385 
386  nPrimaryMatchedPfos.push_back(nPrimaryMatches);
387  nPrimaryMatchedTBHierarchyPfos.push_back(nPrimaryTBHierarchyMatches);
388  nPrimaryMatchedCRPfos.push_back(nPrimaryCRMatches);
389  nTargetMatches += nPrimaryMatches;
390  nTargetTBHierarchyMatches += nPrimaryTBHierarchyMatches;
391  nTargetCRMatches += nPrimaryCRMatches;
392  nTargetGoodTBHierarchyMatches += nPrimaryGoodTBHierarchyMatches;
393  nTargetTBHierarchySplits += nPrimaryTBHierarchySplits;
394  if (0 == nPrimaryMatches)
395  ++nTargetTBHierarchyLosses;
396 
397  if (fillTree)
398  {
399  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "fileIdentifier", m_fileIdentifier));
400  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "eventNumber", eventNumber));
401  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcNuanceCode", mcNuanceCode));
402  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isBeamParticle", isBeamParticle));
403  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isCosmicRay", isCosmicRay));
404  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetPrimaries", nTargetPrimaries));
405  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "targetVertexX", targetVertexX));
406  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "targetVertexY", targetVertexY));
407  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "targetVertexZ", targetVertexZ));
408  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "recoVertexX", recoVertexX));
409  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "recoVertexY", recoVertexY));
410  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "recoVertexZ", recoVertexZ));
411  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryId", &mcPrimaryId));
412  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryPdg", &mcPrimaryPdg));
413  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryTier", &mcPrimaryTier));
414  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryE", &mcPrimaryE));
415  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryPX", &mcPrimaryPX));
416  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryPY", &mcPrimaryPY));
417  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryPZ", &mcPrimaryPZ));
418  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryVtxX", &mcPrimaryVtxX));
419  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryVtxY", &mcPrimaryVtxY));
420  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryVtxZ", &mcPrimaryVtxZ));
421  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryEndX", &mcPrimaryEndX));
422  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryEndY", &mcPrimaryEndY));
423  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryEndZ", &mcPrimaryEndZ));
424  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryNHitsTotal", &nMCHitsTotal));
425  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryNHitsU", &nMCHitsU));
426  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryNHitsV", &nMCHitsV));
427  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryNHitsW", &nMCHitsW));
428  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nPrimaryMatchedPfos", &nPrimaryMatchedPfos));
429  PANDORA_MONITORING_API(
430  SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nPrimaryMatchedTBHierarchyPfos", &nPrimaryMatchedTBHierarchyPfos));
431  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nPrimaryMatchedCRPfos", &nPrimaryMatchedCRPfos));
432  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoId", &bestMatchPfoId));
433  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoPdg", &bestMatchPfoPdg));
434  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoTier", &bestMatchPfoTier));
435  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNHitsTotal", &bestMatchPfoNHitsTotal));
436  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNHitsU", &bestMatchPfoNHitsU));
437  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNHitsV", &bestMatchPfoNHitsV));
438  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNHitsW", &bestMatchPfoNHitsW));
439  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNSharedHitsTotal", &bestMatchPfoNSharedHitsTotal));
440  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNSharedHitsU", &bestMatchPfoNSharedHitsU));
441  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNSharedHitsV", &bestMatchPfoNSharedHitsV));
442  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNSharedHitsW", &bestMatchPfoNSharedHitsW));
443  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoX0", &bestMatchPfoX0));
444  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetMatches", nTargetMatches));
445  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetTBHierarchyMatches", nTargetTBHierarchyMatches));
446  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetCRMatches", nTargetCRMatches));
447 
448  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoIsTestBeam", &bestMatchPfoIsTestBeam));
449  PANDORA_MONITORING_API(
450  SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoIsTestBeamHierarchy", &bestMatchPfoIsTestBeamHierarchy));
451  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoRecoTBId", &bestMatchPfoRecoTBId));
452  PANDORA_MONITORING_API(
453  SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetGoodTBHierarchyMatches", nTargetGoodTBHierarchyMatches));
454  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetTBHierarchySplits", nTargetTBHierarchySplits));
455  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetTBHierarchyLosses", nTargetTBHierarchyLosses));
456  }
457 
458  if (isCosmicRay || isLastTestBeamLeading)
459  {
460  const LArInteractionTypeHelper::InteractionType interactionType(
461  LArInteractionTypeHelper::GetTestBeamHierarchyInteractionType(associatedMCPrimaries));
462 #ifdef MONITORING
463  const int interactionTypeInt(static_cast<int>(interactionType));
464 #endif
465  // ATTN Some redundancy introduced to contributing variables
466  const int isCorrectTBHierarchy(isLeadingBeamParticle && (nTargetGoodTBHierarchyMatches == nTargetTBHierarchyMatches) &&
467  (nTargetGoodTBHierarchyMatches == nTargetPrimaries) && (nTargetCRMatches == 0) && (nTargetTBHierarchySplits == 0) &&
468  (nTargetTBHierarchyLosses == 0));
469  const int isCorrectCR(isCosmicRay && (nTargetTBHierarchyMatches == 0) && (nTargetCRMatches == 1));
470  const int isFakeTBHierarchy(isCosmicRay && (nTargetTBHierarchyMatches > 0));
471  const int isFakeCR(!isCosmicRay && (nTargetCRMatches > 0));
472  const int isSplitTBHierarchy(!isCosmicRay && ((nTargetTBHierarchyMatches > nTargetPrimaries) || (nTargetTBHierarchySplits > 0)));
473  const int isSplitCR(isCosmicRay && (nTargetCRMatches > 1));
474  const int isLost(nTargetMatches == 0);
475 
476  std::stringstream outcomeSS;
477  const bool isBeamHierarchy((mcNuanceCode == 2001) | (mcNuanceCode == 2000));
478  outcomeSS << LArInteractionTypeHelper::ToString(interactionType) << " (Nuance " << mcNuanceCode << ", TB " << isBeamHierarchy
479  << ", CR " << isCosmicRay << ")" << std::endl;
480 
481  if (isLastTestBeamLeading)
482  ++nTotalTBHierarchy;
483  if (isCosmicRay)
484  ++nTotalCR;
485  if (isCorrectTBHierarchy)
486  ++nCorrectTBHierarchy;
487  if (isCorrectCR)
488  ++nCorrectCR;
489  if (isFakeTBHierarchy)
490  ++nFakeTBHierarchy;
491  if (isFakeCR)
492  ++nFakeCR;
493  if (isSplitTBHierarchy)
494  ++nSplitTBHierarchy;
495  if (isSplitCR)
496  ++nSplitCR;
497  if (isLost)
498  ++nLost;
499 
500  if (isCorrectTBHierarchy)
501  outcomeSS << "IsCorrectTBHierarchy";
502  if (isCorrectCR)
503  outcomeSS << "IsCorrectCR ";
504  if (isFakeTBHierarchy)
505  outcomeSS << "IsFake" << name << " ";
506  if (isFakeCR)
507  outcomeSS << "IsFakeCR ";
508  if (isSplitTBHierarchy)
509  outcomeSS << "isSplit" << name << " ";
510  if (isSplitCR)
511  outcomeSS << "IsSplitCR ";
512  if (isLost)
513  outcomeSS << "IsLost ";
514  if (nTargetTBHierarchyMatches > 0)
515  outcomeSS << "(N" << name << "Matches: " << nTargetTBHierarchyMatches << ") ";
516  if (nTargetTBHierarchyLosses > 0)
517  outcomeSS << "(N" << name << "Losses: " << nTargetTBHierarchyLosses << ") ";
518  if (nTargetTBHierarchySplits > 0)
519  outcomeSS << "(N" << name << "Splits: " << nTargetTBHierarchySplits << ") ";
520  if (nTargetCRMatches > 0)
521  outcomeSS << "(NCRMatches: " << nTargetCRMatches << ") ";
522  if (printToScreen)
523  std::cout << outcomeSS.str() << std::endl << targetSS.str() << std::endl;
524 
525  if (fillTree)
526  {
527  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "interactionType", interactionTypeInt));
528  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isCorrectTBHierarchy", isCorrectTBHierarchy));
529  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isCorrectCR", isCorrectCR));
530  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isFakeTBHierarchy", isFakeTBHierarchy));
531  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isFakeCR", isFakeCR));
532  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isSplitTBHierarchy", isSplitTBHierarchy));
533  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isSplitCR", isSplitCR));
534  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isLost", isLost));
535  PANDORA_MONITORING_API(FillTree(this->GetPandora(), m_treeName.c_str()));
536  }
537 
538  targetSS.str(std::string());
539  targetSS.clear();
540  associatedMCPrimaries.clear();
541  nTargetMatches = 0;
542  nTargetTBHierarchyMatches = 0;
543  nTargetCRMatches = 0;
544  nTargetGoodTBHierarchyMatches = 0;
545  nTargetTBHierarchySplits = 0;
546  nTargetTBHierarchyLosses = 0;
547  mcPrimaryId.clear();
548  mcPrimaryPdg.clear();
549  mcPrimaryTier.clear();
550  nMCHitsTotal.clear();
551  nMCHitsU.clear();
552  nMCHitsV.clear();
553  nMCHitsW.clear();
554  mcPrimaryE.clear();
555  mcPrimaryPX.clear();
556  mcPrimaryPY.clear();
557  mcPrimaryPZ.clear();
558  mcPrimaryVtxX.clear();
559  mcPrimaryVtxY.clear();
560  mcPrimaryVtxZ.clear();
561  mcPrimaryEndX.clear();
562  mcPrimaryEndY.clear();
563  mcPrimaryEndZ.clear();
564  nPrimaryMatchedPfos.clear();
565  nPrimaryMatchedTBHierarchyPfos.clear();
566  nPrimaryMatchedCRPfos.clear();
567  bestMatchPfoId.clear();
568  bestMatchPfoPdg.clear();
569  bestMatchPfoTier.clear();
570  bestMatchPfoIsTestBeam.clear();
571  bestMatchPfoIsTestBeamHierarchy.clear();
572  bestMatchPfoRecoTBId.clear();
573  bestMatchPfoNHitsTotal.clear();
574  bestMatchPfoNHitsU.clear();
575  bestMatchPfoNHitsV.clear();
576  bestMatchPfoNHitsW.clear();
577  bestMatchPfoNSharedHitsTotal.clear();
578  bestMatchPfoNSharedHitsU.clear();
579  bestMatchPfoNSharedHitsV.clear();
580  bestMatchPfoNSharedHitsW.clear();
581  bestMatchPfoX0.clear();
582  }
583  }
584 
585  if (useInterpretedMatching)
586  {
587  std::stringstream summarySS;
588  summarySS << "---SUMMARY--------------------------------------------------------------------------------------" << std::endl;
589  if (nTotalTBHierarchy > 0)
590  summarySS << "#CorrectTBHierarchy: " << nCorrectTBHierarchy << "/" << nTotalTBHierarchy << ", Fraction: "
591  << (nTotalTBHierarchy > 0 ? static_cast<float>(nCorrectTBHierarchy) / static_cast<float>(nTotalTBHierarchy) : 0.f)
592  << std::endl;
593  if (nTotalCR > 0)
594  summarySS << "#CorrectCR: " << nCorrectCR << "/" << nTotalCR
595  << ", Fraction: " << (nTotalCR > 0 ? static_cast<float>(nCorrectCR) / static_cast<float>(nTotalCR) : 0.f) << std::endl;
596  if (nFakeTBHierarchy > 0)
597  summarySS << "#Fake" << name << ": " << nFakeTBHierarchy << " ";
598  if (nFakeCR > 0)
599  summarySS << "#FakeCR: " << nFakeCR << " ";
600  if (nSplitTBHierarchy > 0)
601  summarySS << "#Split" << name << ": " << nSplitTBHierarchy << " ";
602  if (nSplitCR > 0)
603  summarySS << "#SplitCR: " << nSplitCR << " ";
604  if (nLost > 0)
605  summarySS << "#Lost: " << nLost << " ";
606  if (nFakeTBHierarchy || nFakeCR || nSplitTBHierarchy || nSplitCR || nLost)
607  summarySS << std::endl;
608  if (printToScreen)
609  std::cout << summarySS.str();
610  }
611 
612  if (printToScreen)
613  std::cout << "------------------------------------------------------------------------------------------------" << std::endl
614  << std::endl;
615 }
616 
617 //------------------------------------------------------------------------------------------------------------------------------------------
618 
619 StatusCode TestBeamHierarchyEventValidationAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
620 {
621  return EventValidationBaseAlgorithm::ReadSettings(xmlHandle);
622 }
623 
624 } // namespace lar_content
const LArMCParticleHelper::PfoContributionMap & GetPfoToHitsMap() const
Get the pfo to hits map.
Header file for the pfo helper class.
std::unordered_map< const pandora::MCParticle *, pandora::CaloHitList > MCContributionMap
const LArMCParticleHelper::MCContributionMap & GetTargetMCParticleToHitsMap() const
Get the target mc particle to hits map.
Header file for the interaction type helper class.
unsigned int m_minPrimaryGoodHits
the minimum number of primary good Hits
const LArMCParticleHelper::MCContributionMap & GetAllMCParticleToHitsMap() const
Get the all mc particle to hits map.
std::map< const pandora::MCParticle *, PfoToSharedHitsVector > MCParticleToPfoHitSharingMap
const LArMCParticleHelper::MCParticleToPfoHitSharingMap & GetMCToPfoHitSharingMap() const
Get the mc to pfo hit sharing map.
void SetInterpretedMCToPfoHitSharingMap(const LArMCParticleHelper::MCParticleToPfoHitSharingMap &interpretedMCToPfoHitSharingMap)
Set the interpreted mc to pfo hit sharing map.
std::vector< int > IntVector
Header file for the lar monitoring helper helper class.
void SetMCToPfoHitSharingMap(const LArMCParticleHelper::MCParticleToPfoHitSharingMap &mcToPfoHitSharingMap)
Set the mc to pfo hit sharing map.
unsigned int m_minHitsForGoodView
the minimum number of Hits for a good view
Header file for the test beam hierarchy event validation algorithm.
void SetTargetMCParticleToHitsMap(const LArMCParticleHelper::MCContributionMap &targetMCParticleToHitsMap)
Set the target mc particle to hits map.
std::unordered_map< const pandora::MCParticle *, int > MCParticleIntMap
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
const LArMCParticleHelper::MCParticleToPfoHitSharingMap & GetInterpretedMCToPfoHitSharingMap() const
Get the interpreted mc to pfo hit sharing map.
float m_minHitSharingFraction
the minimum Hit sharing fraction
void SetAllMCParticleToHitsMap(const LArMCParticleHelper::MCContributionMap &allMCParticleToHitsMap)
Set the all mc particle to hits map.
std::pair< const pandora::ParticleFlowObject *, pandora::CaloHitList > PfoCaloHitListPair
void SetPfoToHitsMap(const LArMCParticleHelper::PfoContributionMap &pfoToHitsMap)
Set the pfo to hits map.
std::unordered_map< const pandora::ParticleFlowObject *, unsigned int > PfoToIdMap
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
std::unordered_map< const pandora::ParticleFlowObject *, pandora::CaloHitList > PfoContributionMap
std::unordered_map< const pandora::MCParticle *, const pandora::MCParticle * > MCRelationMap
std::map< const pandora::ParticleFlowObject *, MCParticleToSharedHitsVector > PfoToMCParticleHitSharingMap