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