All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
NeutrinoEventValidationAlgorithm.cc
Go to the documentation of this file.
1 /**
2  * @file larpandoracontent/LArMonitoring/NeutrinoEventValidationAlgorithm.cc
3  *
4  * @brief Implementation of the neutrino event validation algorithm.
5  *
6  * $Log: $
7  */
8 
9 #include "Pandora/AlgorithmHeaders.h"
10 
14 
16 
17 #include <sstream>
18 
19 using namespace pandora;
20 
21 namespace lar_content
22 {
23 
24 NeutrinoEventValidationAlgorithm::NeutrinoEventValidationAlgorithm() : m_useTrueNeutrinosOnly(false)
25 {
26 }
27 
28 //------------------------------------------------------------------------------------------------------------------------------------------
29 
31 {
32 }
33 
34 //------------------------------------------------------------------------------------------------------------------------------------------
35 
36 void NeutrinoEventValidationAlgorithm::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;
43  pMCParticleList, pCaloHitList, m_primaryParameters, LArMCParticleHelper::IsBeamNeutrinoFinalState, targetMCParticleToHitsMap);
46  pMCParticleList, pCaloHitList, m_primaryParameters, LArMCParticleHelper::IsCosmicRay, targetMCParticleToHitsMap);
47 
49  parameters.m_minPrimaryGoodHits = 0;
50  parameters.m_minHitsForGoodView = 0;
51  parameters.m_minHitSharingFraction = 0.f;
52  LArMCParticleHelper::MCContributionMap allMCParticleToHitsMap;
54  pMCParticleList, pCaloHitList, parameters, LArMCParticleHelper::IsBeamNeutrinoFinalState, allMCParticleToHitsMap);
57  pMCParticleList, pCaloHitList, parameters, LArMCParticleHelper::IsCosmicRay, allMCParticleToHitsMap);
58 
59  validationInfo.SetTargetMCParticleToHitsMap(targetMCParticleToHitsMap);
60  validationInfo.SetAllMCParticleToHitsMap(allMCParticleToHitsMap);
61  }
62 
63  if (pPfoList)
64  {
65  PfoList allConnectedPfos;
66  LArPfoHelper::GetAllConnectedPfos(*pPfoList, allConnectedPfos);
67 
68  PfoList finalStatePfos;
69  for (const ParticleFlowObject *const pPfo : allConnectedPfos)
70  {
72  finalStatePfos.push_back(pPfo);
73  }
74 
77  finalStatePfos, validationInfo.GetAllMCParticleToHitsMap(), pfoToHitsMap, m_primaryParameters.m_foldBackHierarchy);
78 
79  validationInfo.SetPfoToHitsMap(pfoToHitsMap);
80  }
81 
85  validationInfo.GetPfoToHitsMap(), {validationInfo.GetAllMCParticleToHitsMap()}, pfoToMCHitSharingMap, mcToPfoHitSharingMap);
86 
87  // ATTN : Ensure all mc primaries have an entry in mcToPfoHitSharingMap, even if no associated pfos.
88  MCParticleVector mcPrimaryVector;
90  for (const MCParticle *pMCParticle : mcPrimaryVector)
91  {
92  if (mcToPfoHitSharingMap.find(pMCParticle) == mcToPfoHitSharingMap.end())
93  {
94  LArMCParticleHelper::PfoToSharedHitsVector pfoToSharedHitsVector;
95  mcToPfoHitSharingMap.insert(LArMCParticleHelper::MCParticleToPfoHitSharingMap::value_type(pMCParticle, pfoToSharedHitsVector));
96  }
97  }
98 
99  validationInfo.SetMCToPfoHitSharingMap(mcToPfoHitSharingMap);
100 
101  LArMCParticleHelper::MCParticleToPfoHitSharingMap interpretedMCToPfoHitSharingMap;
102  this->InterpretMatching(validationInfo, interpretedMCToPfoHitSharingMap);
103  validationInfo.SetInterpretedMCToPfoHitSharingMap(interpretedMCToPfoHitSharingMap);
104 }
105 
106 //------------------------------------------------------------------------------------------------------------------------------------------
107 
109  const ValidationInfo &validationInfo, const bool useInterpretedMatching, const bool printToScreen, const bool fillTree) const
110 {
111  if (printToScreen && useInterpretedMatching)
112  std::cout << "---INTERPRETED-MATCHING-OUTPUT------------------------------------------------------------------" << std::endl;
113  else if (printToScreen)
114  std::cout << "---RAW-MATCHING-OUTPUT--------------------------------------------------------------------------" << std::endl;
115 
116  const LArMCParticleHelper::MCParticleToPfoHitSharingMap &mcToPfoHitSharingMap(
117  useInterpretedMatching ? validationInfo.GetInterpretedMCToPfoHitSharingMap() : validationInfo.GetMCToPfoHitSharingMap());
118 
119  MCParticleVector mcPrimaryVector;
121 
122  // Neutrino Validation Bookkeeping
123  int nNeutrinoPrimaries(0);
124  for (const MCParticle *const pMCPrimary : mcPrimaryVector)
125  if (LArMCParticleHelper::IsBeamNeutrinoFinalState(pMCPrimary) && validationInfo.GetTargetMCParticleToHitsMap().count(pMCPrimary))
126  ++nNeutrinoPrimaries;
127 
128  PfoVector primaryPfoVector;
129  LArMonitoringHelper::GetOrderedPfoVector(validationInfo.GetPfoToHitsMap(), primaryPfoVector);
130 
131  int pfoIndex(0), neutrinoPfoIndex(0);
132  PfoToIdMap pfoToIdMap, neutrinoPfoToIdMap;
133 
134  for (const Pfo *const pPrimaryPfo : primaryPfoVector)
135  {
136  pfoToIdMap.insert(PfoToIdMap::value_type(pPrimaryPfo, ++pfoIndex));
137  const Pfo *const pRecoNeutrino(LArPfoHelper::IsNeutrinoFinalState(pPrimaryPfo) ? LArPfoHelper::GetParentNeutrino(pPrimaryPfo) : nullptr);
138 
139  if (pRecoNeutrino && !neutrinoPfoToIdMap.count(pRecoNeutrino))
140  neutrinoPfoToIdMap.insert(PfoToIdMap::value_type(pRecoNeutrino, ++neutrinoPfoIndex));
141  }
142 
143  LArMCParticleHelper::MCParticleIntMap triggeredToLeading, triggeredToLeadingCounter;
144 
145  PfoSet recoNeutrinos;
146  MCParticleList associatedMCPrimaries;
147 
148  int nCorrectNu(0), nTotalNu(0), nCorrectCR(0), nTotalCR(0);
149  int nFakeNu(0), nFakeCR(0), nSplitNu(0), nSplitCR(0), nLost(0), mcPrimaryIndex(0), nTargetMatches(0), nTargetNuMatches(0);
150  int nTargetCRMatches(0), nTargetGoodNuMatches(0), nTargetNuSplits(0), nTargetNuLosses(0);
151  IntVector mcPrimaryId, mcPrimaryPdg, nMCHitsTotal, nMCHitsU, nMCHitsV, nMCHitsW;
152  FloatVector mcPrimaryE, mcPrimaryPX, mcPrimaryPY, mcPrimaryPZ;
153  FloatVector mcPrimaryVtxX, mcPrimaryVtxY, mcPrimaryVtxZ, mcPrimaryEndX, mcPrimaryEndY, mcPrimaryEndZ;
154  IntVector nPrimaryMatchedPfos, nPrimaryMatchedNuPfos, nPrimaryMatchedCRPfos;
155  IntVector bestMatchPfoId, bestMatchPfoPdg, bestMatchPfoIsRecoNu, bestMatchPfoRecoNuId;
156  IntVector bestMatchPfoNHitsTotal, bestMatchPfoNHitsU, bestMatchPfoNHitsV, bestMatchPfoNHitsW;
157  IntVector bestMatchPfoNSharedHitsTotal, bestMatchPfoNSharedHitsU, bestMatchPfoNSharedHitsV, bestMatchPfoNSharedHitsW;
158 
159  std::stringstream targetSS;
160  const std::string name("Nu");
161 
162  for (const MCParticle *const pMCPrimary : mcPrimaryVector)
163  {
164  const bool hasMatch(mcToPfoHitSharingMap.count(pMCPrimary) && !mcToPfoHitSharingMap.at(pMCPrimary).empty());
165  const bool isTargetPrimary(validationInfo.GetTargetMCParticleToHitsMap().count(pMCPrimary));
166 
167  if (!isTargetPrimary && !hasMatch)
168  continue;
169 
170  associatedMCPrimaries.push_back(pMCPrimary);
171  const int nTargetPrimaries(associatedMCPrimaries.size());
172  const bool isLastNeutrinoPrimary(++mcPrimaryIndex == nNeutrinoPrimaries);
173  const CaloHitList &mcPrimaryHitList(validationInfo.GetAllMCParticleToHitsMap().at(pMCPrimary));
174 
176  const int isBeamNeutrinoFinalState(LArMCParticleHelper::IsBeamNeutrinoFinalState(pMCPrimary));
177  const int isCosmicRay(LArMCParticleHelper::IsCosmicRay(pMCPrimary));
178 #ifdef MONITORING
179  const CartesianVector &targetVertex(LArMCParticleHelper::GetParentMCParticle(pMCPrimary)->GetVertex());
180  const float targetVertexX(targetVertex.GetX()), targetVertexY(targetVertex.GetY()), targetVertexZ(targetVertex.GetZ());
181 #endif
182 
183  targetSS << (!isTargetPrimary ? "(Non target) " : "") << "PrimaryId " << mcPrimaryIndex << ", Nu " << isBeamNeutrinoFinalState
184  << ", CR " << isCosmicRay << ", MCPDG " << pMCPrimary->GetParticleId() << ", Energy " << pMCPrimary->GetEnergy()
185  << ", Dist. " << (pMCPrimary->GetEndpoint() - pMCPrimary->GetVertex()).GetMagnitude() << ", nMCHits "
186  << mcPrimaryHitList.size() << " (" << LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, mcPrimaryHitList) << ", "
187  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, mcPrimaryHitList) << ", "
188  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, mcPrimaryHitList) << ")" << std::endl;
189 
190  mcPrimaryId.push_back(mcPrimaryIndex);
191  mcPrimaryPdg.push_back(pMCPrimary->GetParticleId());
192  mcPrimaryE.push_back(pMCPrimary->GetEnergy());
193  mcPrimaryPX.push_back(pMCPrimary->GetMomentum().GetX());
194  mcPrimaryPY.push_back(pMCPrimary->GetMomentum().GetY());
195  mcPrimaryPZ.push_back(pMCPrimary->GetMomentum().GetZ());
196  mcPrimaryVtxX.push_back(pMCPrimary->GetVertex().GetX());
197  mcPrimaryVtxY.push_back(pMCPrimary->GetVertex().GetY());
198  mcPrimaryVtxZ.push_back(pMCPrimary->GetVertex().GetZ());
199  mcPrimaryEndX.push_back(pMCPrimary->GetEndpoint().GetX());
200  mcPrimaryEndY.push_back(pMCPrimary->GetEndpoint().GetY());
201  mcPrimaryEndZ.push_back(pMCPrimary->GetEndpoint().GetZ());
202  nMCHitsTotal.push_back(mcPrimaryHitList.size());
203  nMCHitsU.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, mcPrimaryHitList));
204  nMCHitsV.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, mcPrimaryHitList));
205  nMCHitsW.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, mcPrimaryHitList));
206 
207  int matchIndex(0), nPrimaryMatches(0), nPrimaryNuMatches(0), nPrimaryCRMatches(0), nPrimaryGoodNuMatches(0), nPrimaryNuSplits(0);
208 #ifdef MONITORING
209  float recoVertexX(std::numeric_limits<float>::max()), recoVertexY(std::numeric_limits<float>::max()),
210  recoVertexZ(std::numeric_limits<float>::max());
211 #endif
212  for (const LArMCParticleHelper::PfoCaloHitListPair &pfoToSharedHits : mcToPfoHitSharingMap.at(pMCPrimary))
213  {
214  const CaloHitList &sharedHitList(pfoToSharedHits.second);
215  const CaloHitList &pfoHitList(validationInfo.GetPfoToHitsMap().at(pfoToSharedHits.first));
216 
217  const bool isRecoNeutrinoFinalState(LArPfoHelper::IsNeutrinoFinalState(pfoToSharedHits.first));
218  const bool isGoodMatch(this->IsGoodMatch(mcPrimaryHitList, pfoHitList, sharedHitList));
219 
220  const int pfoId(pfoToIdMap.at(pfoToSharedHits.first));
221  const int recoNuId(isRecoNeutrinoFinalState ? neutrinoPfoToIdMap.at(LArPfoHelper::GetParentNeutrino(pfoToSharedHits.first)) : -1);
222 
223  if (0 == matchIndex++)
224  {
225  bestMatchPfoId.push_back(pfoId);
226  bestMatchPfoPdg.push_back(pfoToSharedHits.first->GetParticleId());
227  bestMatchPfoIsRecoNu.push_back(isRecoNeutrinoFinalState ? 1 : 0);
228  bestMatchPfoRecoNuId.push_back(recoNuId);
229  bestMatchPfoNHitsTotal.push_back(pfoHitList.size());
230  bestMatchPfoNHitsU.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, pfoHitList));
231  bestMatchPfoNHitsV.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, pfoHitList));
232  bestMatchPfoNHitsW.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, pfoHitList));
233  bestMatchPfoNSharedHitsTotal.push_back(sharedHitList.size());
234  bestMatchPfoNSharedHitsU.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, sharedHitList));
235  bestMatchPfoNSharedHitsV.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, sharedHitList));
236  bestMatchPfoNSharedHitsW.push_back(LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, sharedHitList));
237 #ifdef MONITORING
238  try
239  {
240  const Vertex *const pRecoVertex(LArPfoHelper::GetVertex(
241  isRecoNeutrinoFinalState ? LArPfoHelper::GetParentNeutrino(pfoToSharedHits.first) : pfoToSharedHits.first));
242  recoVertexX = pRecoVertex->GetPosition().GetX();
243  recoVertexY = pRecoVertex->GetPosition().GetY();
244  recoVertexZ = pRecoVertex->GetPosition().GetZ();
245  }
246  catch (const StatusCodeException &)
247  {
248  }
249 #endif
250  }
251 
252  if (isGoodMatch)
253  ++nPrimaryMatches;
254 
255  if (isRecoNeutrinoFinalState)
256  {
257  const Pfo *const pRecoNeutrino(LArPfoHelper::GetParentNeutrino(pfoToSharedHits.first));
258  const bool isSplitRecoNeutrino(!recoNeutrinos.empty() && !recoNeutrinos.count(pRecoNeutrino));
259  if (!isSplitRecoNeutrino && isGoodMatch)
260  ++nPrimaryGoodNuMatches;
261  if (isSplitRecoNeutrino && isBeamNeutrinoFinalState && isGoodMatch)
262  ++nPrimaryNuSplits;
263  recoNeutrinos.insert(pRecoNeutrino);
264  }
265 
266  if (isRecoNeutrinoFinalState && isGoodMatch)
267  ++nPrimaryNuMatches;
268  if (!isRecoNeutrinoFinalState && isGoodMatch)
269  ++nPrimaryCRMatches;
270 
271  targetSS << "-" << (!isGoodMatch ? "(Below threshold) " : "") << "MatchedPfoId " << pfoId << ", Nu " << isRecoNeutrinoFinalState;
272  if (isRecoNeutrinoFinalState)
273  targetSS << " [NuId: " << recoNuId << "]";
274  targetSS << ", CR " << (!isRecoNeutrinoFinalState) << ", PDG " << pfoToSharedHits.first->GetParticleId() << ", nMatchedHits "
275  << sharedHitList.size() << " (" << LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, sharedHitList) << ", "
276  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, sharedHitList) << ", "
277  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, sharedHitList) << ")"
278  << ", nPfoHits " << pfoHitList.size() << " (" << LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, pfoHitList) << ", "
279  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, pfoHitList) << ", "
280  << LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, pfoHitList) << ")" << std::endl;
281  }
282 
283  if (mcToPfoHitSharingMap.at(pMCPrimary).empty())
284  {
285  targetSS << "-No matched Pfo" << std::endl;
286  bestMatchPfoId.push_back(-1);
287  bestMatchPfoPdg.push_back(0);
288  bestMatchPfoIsRecoNu.push_back(0);
289  bestMatchPfoRecoNuId.push_back(-1);
290  bestMatchPfoNHitsTotal.push_back(0);
291  bestMatchPfoNHitsU.push_back(0);
292  bestMatchPfoNHitsV.push_back(0);
293  bestMatchPfoNHitsW.push_back(0);
294  bestMatchPfoNSharedHitsTotal.push_back(0);
295  bestMatchPfoNSharedHitsU.push_back(0);
296  bestMatchPfoNSharedHitsV.push_back(0);
297  bestMatchPfoNSharedHitsW.push_back(0);
298  }
299 
300  nPrimaryMatchedPfos.push_back(nPrimaryMatches);
301  nPrimaryMatchedNuPfos.push_back(nPrimaryNuMatches);
302  nPrimaryMatchedCRPfos.push_back(nPrimaryCRMatches);
303  nTargetMatches += nPrimaryMatches;
304  nTargetNuMatches += nPrimaryNuMatches;
305  nTargetCRMatches += nPrimaryCRMatches;
306  nTargetGoodNuMatches += nPrimaryGoodNuMatches;
307  nTargetNuSplits += nPrimaryNuSplits;
308  if (0 == nPrimaryMatches)
309  ++nTargetNuLosses;
310 
311  if (fillTree)
312  {
313  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "fileIdentifier", m_fileIdentifier));
314  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "eventNumber", m_eventNumber - 1));
315  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcNuanceCode", mcNuanceCode));
316  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isNeutrino", isBeamNeutrinoFinalState));
317  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isCosmicRay", isCosmicRay));
318  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetPrimaries", nTargetPrimaries));
319  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "targetVertexX", targetVertexX));
320  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "targetVertexY", targetVertexY));
321  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "targetVertexZ", targetVertexZ));
322  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "recoVertexX", recoVertexX));
323  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "recoVertexY", recoVertexY));
324  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "recoVertexZ", recoVertexZ));
325  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryId", &mcPrimaryId));
326  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryPdg", &mcPrimaryPdg));
327  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryE", &mcPrimaryE));
328  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryPX", &mcPrimaryPX));
329  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryPY", &mcPrimaryPY));
330  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryPZ", &mcPrimaryPZ));
331  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryVtxX", &mcPrimaryVtxX));
332  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryVtxY", &mcPrimaryVtxY));
333  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryVtxZ", &mcPrimaryVtxZ));
334  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryEndX", &mcPrimaryEndX));
335  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryEndY", &mcPrimaryEndY));
336  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryEndZ", &mcPrimaryEndZ));
337  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryNHitsTotal", &nMCHitsTotal));
338  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryNHitsU", &nMCHitsU));
339  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryNHitsV", &nMCHitsV));
340  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "mcPrimaryNHitsW", &nMCHitsW));
341  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nPrimaryMatchedPfos", &nPrimaryMatchedPfos));
342  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nPrimaryMatchedNuPfos", &nPrimaryMatchedNuPfos));
343  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nPrimaryMatchedCRPfos", &nPrimaryMatchedCRPfos));
344  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoId", &bestMatchPfoId));
345  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoPdg", &bestMatchPfoPdg));
346  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoIsRecoNu", &bestMatchPfoIsRecoNu));
347  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoRecoNuId", &bestMatchPfoRecoNuId));
348  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNHitsTotal", &bestMatchPfoNHitsTotal));
349  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNHitsU", &bestMatchPfoNHitsU));
350  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNHitsV", &bestMatchPfoNHitsV));
351  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNHitsW", &bestMatchPfoNHitsW));
352  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNSharedHitsTotal", &bestMatchPfoNSharedHitsTotal));
353  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNSharedHitsU", &bestMatchPfoNSharedHitsU));
354  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNSharedHitsV", &bestMatchPfoNSharedHitsV));
355  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "bestMatchPfoNSharedHitsW", &bestMatchPfoNSharedHitsW));
356  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetMatches", nTargetMatches));
357  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetNuMatches", nTargetNuMatches));
358  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetCRMatches", nTargetCRMatches));
359  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetGoodNuMatches", nTargetGoodNuMatches));
360  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetNuSplits", nTargetNuSplits));
361  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "nTargetNuLosses", nTargetNuLosses));
362  }
363 
364  if (isLastNeutrinoPrimary || isCosmicRay)
365  {
367 #ifdef MONITORING
368  const int interactionTypeInt(static_cast<int>(interactionType));
369 #endif
370  // ATTN Some redundancy introduced to contributing variables
371  const int isCorrectNu(isBeamNeutrinoFinalState && (nTargetGoodNuMatches == nTargetNuMatches) && (nTargetGoodNuMatches == nTargetPrimaries) &&
372  (nTargetCRMatches == 0) && (nTargetNuSplits == 0) && (nTargetNuLosses == 0));
373  const int isCorrectCR(isCosmicRay && (nTargetNuMatches == 0) && (nTargetCRMatches == 1));
374  const int isFakeNu(isCosmicRay && (nTargetNuMatches > 0));
375  const int isFakeCR(!isCosmicRay && (nTargetCRMatches > 0));
376  const int isSplitNu(!isCosmicRay && ((nTargetNuMatches > nTargetPrimaries) || (nTargetNuSplits > 0)));
377  const int isSplitCR(isCosmicRay && (nTargetCRMatches > 1));
378  const int isLost(nTargetMatches == 0);
379 
380  std::stringstream outcomeSS;
381  outcomeSS << LArInteractionTypeHelper::ToString(interactionType) << " (Nuance " << mcNuanceCode << ", Nu "
382  << isBeamNeutrinoFinalState << ", CR " << isCosmicRay << ")" << std::endl;
383 
384  if (isLastNeutrinoPrimary)
385  ++nTotalNu;
386  if (isCosmicRay)
387  ++nTotalCR;
388  if (isCorrectNu)
389  ++nCorrectNu;
390  if (isCorrectCR)
391  ++nCorrectCR;
392  if (isFakeNu)
393  ++nFakeNu;
394  if (isFakeCR)
395  ++nFakeCR;
396  if (isSplitNu)
397  ++nSplitNu;
398  if (isSplitCR)
399  ++nSplitCR;
400  if (isLost)
401  ++nLost;
402 
403  if (isCorrectNu)
404  outcomeSS << "IsCorrectNu ";
405  if (isCorrectCR)
406  outcomeSS << "IsCorrectCR ";
407  if (isFakeNu)
408  outcomeSS << "IsFake" << name << " ";
409  if (isFakeCR)
410  outcomeSS << "IsFakeCR ";
411  if (isSplitNu)
412  outcomeSS << "isSplit" << name << " ";
413  if (isSplitCR)
414  outcomeSS << "IsSplitCR ";
415  if (isLost)
416  outcomeSS << "IsLost ";
417  if (nTargetNuMatches > 0)
418  outcomeSS << "(N" << name << "Matches: " << nTargetNuMatches << ") ";
419  if (nTargetNuLosses > 0)
420  outcomeSS << "(N" << name << "Losses: " << nTargetNuLosses << ") ";
421  if (nTargetNuSplits > 0)
422  outcomeSS << "(N" << name << "Splits: " << nTargetNuSplits << ") ";
423  if (nTargetCRMatches > 0)
424  outcomeSS << "(NCRMatches: " << nTargetCRMatches << ") ";
425  if (printToScreen)
426  std::cout << outcomeSS.str() << std::endl << targetSS.str() << std::endl;
427 
428  if (fillTree)
429  {
430  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "interactionType", interactionTypeInt));
431  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isCorrectNu", isCorrectNu));
432  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isCorrectCR", isCorrectCR));
433  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isFakeNu", isFakeNu));
434  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isFakeCR", isFakeCR));
435  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isSplitNu", isSplitNu));
436  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isSplitCR", isSplitCR));
437  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treeName.c_str(), "isLost", isLost));
438  PANDORA_MONITORING_API(FillTree(this->GetPandora(), m_treeName.c_str()));
439  }
440 
441  targetSS.str(std::string());
442  targetSS.clear();
443  recoNeutrinos.clear();
444  associatedMCPrimaries.clear();
445  nTargetMatches = 0;
446  nTargetNuMatches = 0;
447  nTargetCRMatches = 0;
448  nTargetGoodNuMatches = 0;
449  nTargetNuSplits = 0;
450  nTargetNuLosses = 0;
451  mcPrimaryId.clear();
452  mcPrimaryPdg.clear();
453  nMCHitsTotal.clear();
454  nMCHitsU.clear();
455  nMCHitsV.clear();
456  nMCHitsW.clear();
457  mcPrimaryE.clear();
458  mcPrimaryPX.clear();
459  mcPrimaryPY.clear();
460  mcPrimaryPZ.clear();
461  mcPrimaryVtxX.clear();
462  mcPrimaryVtxY.clear();
463  mcPrimaryVtxZ.clear();
464  mcPrimaryEndX.clear();
465  mcPrimaryEndY.clear();
466  mcPrimaryEndZ.clear();
467  nPrimaryMatchedPfos.clear();
468  nPrimaryMatchedNuPfos.clear();
469  nPrimaryMatchedCRPfos.clear();
470  bestMatchPfoId.clear();
471  bestMatchPfoPdg.clear();
472  bestMatchPfoIsRecoNu.clear();
473  bestMatchPfoRecoNuId.clear();
474  bestMatchPfoNHitsTotal.clear();
475  bestMatchPfoNHitsU.clear();
476  bestMatchPfoNHitsV.clear();
477  bestMatchPfoNHitsW.clear();
478  bestMatchPfoNSharedHitsTotal.clear();
479  bestMatchPfoNSharedHitsU.clear();
480  bestMatchPfoNSharedHitsV.clear();
481  bestMatchPfoNSharedHitsW.clear();
482  }
483  }
484 
485  if (useInterpretedMatching)
486  {
487  std::stringstream summarySS;
488  summarySS << "---SUMMARY--------------------------------------------------------------------------------------" << std::endl;
489  if (nTotalNu > 0)
490  summarySS << "#CorrectNu: " << nCorrectNu << "/" << nTotalNu
491  << ", Fraction: " << (nTotalNu > 0 ? static_cast<float>(nCorrectNu) / static_cast<float>(nTotalNu) : 0.f) << std::endl;
492  if (nTotalCR > 0)
493  summarySS << "#CorrectCR: " << nCorrectCR << "/" << nTotalCR
494  << ", Fraction: " << (nTotalCR > 0 ? static_cast<float>(nCorrectCR) / static_cast<float>(nTotalCR) : 0.f) << std::endl;
495  if (nFakeNu > 0)
496  summarySS << "#Fake" << name << ": " << nFakeNu << " ";
497  if (nFakeCR > 0)
498  summarySS << "#FakeCR: " << nFakeCR << " ";
499  if (nSplitNu > 0)
500  summarySS << "#Split" << name << ": " << nSplitNu << " ";
501  if (nSplitCR > 0)
502  summarySS << "#SplitCR: " << nSplitCR << " ";
503  if (nLost > 0)
504  summarySS << "#Lost: " << nLost << " ";
505  if (nFakeNu || nFakeCR || nSplitNu || nSplitCR || nLost)
506  summarySS << std::endl;
507  if (printToScreen)
508  std::cout << summarySS.str();
509  }
510 
511  if (printToScreen)
512  std::cout << "------------------------------------------------------------------------------------------------" << std::endl
513  << std::endl;
514 }
515 
516 //------------------------------------------------------------------------------------------------------------------------------------------
517 
518 StatusCode NeutrinoEventValidationAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
519 {
520  PANDORA_RETURN_RESULT_IF_AND_IF(
521  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "UseTrueNeutrinosOnly", m_useTrueNeutrinosOnly));
522 
524 }
525 
526 } // namespace lar_content
const LArMCParticleHelper::PfoContributionMap & GetPfoToHitsMap() const
Get the pfo to hits map.
Header file for the pfo helper class.
Header file for the neutrino event validation algorithm.
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
static const pandora::Vertex * GetVertex(const pandora::ParticleFlowObject *const pPfo)
Get the pfo vertex.
void InterpretMatching(const ValidationInfo &validationInfo, LArMCParticleHelper::MCParticleToPfoHitSharingMap &interpretedMCToPfoHitSharingMap) const
Apply an interpretative matching procedure to the comprehensive matches in the provided validation in...
static void GetPfoToReconstructable2DHitsMap(const pandora::PfoList &pfoList, const MCContributionMap &selectedMCParticleToHitsMap, PfoContributionMap &pfoToReconstructable2DHitsMap, const bool foldBackHierarchy)
Get mapping from Pfo to reconstructable 2D hits (=good hits belonging to a selected reconstructable M...
const LArMCParticleHelper::MCContributionMap & GetAllMCParticleToHitsMap() const
Get the all mc particle to hits map.
static void GetPfoMCParticleHitSharingMaps(const PfoContributionMap &pfoToReconstructable2DHitsMap, const MCContributionMapVector &selectedMCParticleToHitsMaps, PfoToMCParticleHitSharingMap &pfoToMCParticleHitSharingMap, MCParticleToPfoHitSharingMap &mcParticleToPfoHitSharingMap)
Get the mappings from Pfo -&gt; pair (reconstructable MCparticles, number of reconstructable 2D hits sha...
std::map< const pandora::MCParticle *, PfoToSharedHitsVector > MCParticleToPfoHitSharingMap
void FillValidationInfo(const pandora::MCParticleList *const pMCParticleList, const pandora::CaloHitList *const pCaloHitList, const pandora::PfoList *const pPfoList, ValidationInfo &validationInfo) const
Fill the validation info containers.
std::vector< PfoCaloHitListPair > PfoToSharedHitsVector
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
bool m_foldBackHierarchy
whether to fold the hierarchy back to the primary (neutrino) or leading particles (test beam) ...
bool m_useTrueNeutrinosOnly
Whether to consider only mc particles that were neutrino induced.
Header file for the lar monitoring helper helper class.
void SetMCToPfoHitSharingMap(const LArMCParticleHelper::MCParticleToPfoHitSharingMap &mcToPfoHitSharingMap)
Set the mc to pfo hit sharing map.
static unsigned int GetNuanceCode(const pandora::MCParticle *const pMCParticle)
Get the nuance code of an MCParticle.
static InteractionType GetInteractionType(const pandora::MCParticleList &mcPrimaryList)
Get the interaction type of an event.
unsigned int m_minHitsForGoodView
the minimum number of Hits for a good view
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
static void GetOrderedMCParticleVector(const LArMCParticleHelper::MCContributionMapVector &selectedMCParticleToGoodHitsMaps, pandora::MCParticleVector &orderedMCParticleVector)
Order input MCParticles by their number of hits.
static bool IsCosmicRay(const pandora::MCParticle *const pMCParticle)
Return true if passed a primary cosmic ray MCParticle.
void SetTargetMCParticleToHitsMap(const LArMCParticleHelper::MCContributionMap &targetMCParticleToHitsMap)
Set the target mc particle to hits map.
static void SelectReconstructableMCParticles(const pandora::MCParticleList *pMCParticleList, const pandora::CaloHitList *pCaloHitList, const PrimaryParameters &parameters, std::function< bool(const pandora::MCParticle *const)> fCriteria, MCContributionMap &selectedMCParticlesToHitsMap)
Select target, reconstructable mc particles that match given criteria.
std::unordered_map< const pandora::MCParticle *, int > MCParticleIntMap
LArMCParticleHelper::PrimaryParameters m_primaryParameters
The mc particle primary selection parameters.
static std::string ToString(const InteractionType interactionType)
Get a string representation of an interaction type.
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
static const pandora::MCParticle * GetParentMCParticle(const pandora::MCParticle *const pMCParticle)
Get the parent mc particle.
static void GetOrderedPfoVector(const LArMCParticleHelper::PfoContributionMap &pfoToReconstructable2DHitsMap, pandora::PfoVector &orderedPfoVector)
Order input Pfos by their number of hits.
static const pandora::ParticleFlowObject * GetParentNeutrino(const pandora::ParticleFlowObject *const pPfo)
Get primary neutrino or antineutrino.
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::unordered_map< const pandora::ParticleFlowObject *, unsigned int > PfoToIdMap
required by fuzzyCluster table::sbnd_g4_services gaushitTruthMatch pandora
Definition: reco_sbnd.fcl:182
std::pair< const pandora::ParticleFlowObject *, pandora::CaloHitList > PfoCaloHitListPair
static bool IsFinalState(const pandora::ParticleFlowObject *const pPfo)
Whether a pfo is a primary parent particle.
void SetPfoToHitsMap(const LArMCParticleHelper::PfoContributionMap &pfoToHitsMap)
Set the pfo to hits map.
static bool IsNeutrinoFinalState(const pandora::ParticleFlowObject *const pPfo)
Whether a pfo is a final-state particle from a neutrino (or antineutrino) interaction.
then echo fcl name
bool IsGoodMatch(const pandora::CaloHitList &trueHits, const pandora::CaloHitList &recoHits, const pandora::CaloHitList &sharedHits) const
Whether a provided mc primary and pfo are deemed to be a good match.
std::unordered_map< const pandora::ParticleFlowObject *, pandora::CaloHitList > PfoContributionMap
static void GetAllConnectedPfos(const pandora::PfoList &inputPfoList, pandora::PfoList &outputPfoList)
Get a flat list of all pfos, recursively including all daughters and parents associated with those pf...
static bool IsBeamNeutrinoFinalState(const pandora::MCParticle *const pMCParticle)
Returns true if passed a primary neutrino final state MCParticle.
static unsigned int CountHitsByType(const pandora::HitType hitType, const pandora::CaloHitList &caloHitList)
Count the number of calo hits, in a provided list, of a specified type.
void ProcessOutput(const ValidationInfo &validationInfo, const bool useInterpretedMatching, const bool printToScreen, const bool fillTree) const
Print matching information in a provided validation info object, and write information to tree if con...
BEGIN_PROLOG could also be cout
std::map< const pandora::ParticleFlowObject *, MCParticleToSharedHitsVector > PfoToMCParticleHitSharingMap
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)