9 #include "Pandora/AlgorithmHeaders.h" 24 EventValidationAlgorithm::EventValidationAlgorithm() :
25 m_useTrueNeutrinosOnly(false),
26 m_testBeamMode(false),
27 m_selectInputHits(true),
28 m_minHitSharingFraction(0.9
f),
29 m_maxPhotonPropagation(2.5
f),
30 m_printAllToScreen(false),
31 m_printMatchingToScreen(true),
33 m_useSmallPrimaries(true),
34 m_matchingMinSharedHits(5),
35 m_matchingMinCompleteness(0.1
f),
36 m_matchingMinPurity(0.5
f),
50 PANDORA_MONITORING_API(SaveTree(this->GetPandora(),
m_treeName.c_str(),
m_fileName.c_str(),
"UPDATE"));
52 catch (
const StatusCodeException &)
54 std::cout <<
"EventValidationAlgorithm: Unable to write tree " <<
m_treeName <<
" to file " <<
m_fileName << std::endl;
65 const MCParticleList *pMCParticleList =
nullptr;
66 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*
this,
m_mcParticleListName, pMCParticleList));
68 const CaloHitList *pCaloHitList =
nullptr;
69 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*
this,
m_caloHitListName, pCaloHitList));
71 const PfoList *pPfoList =
nullptr;
72 (void) PandoraContentApi::GetList(*
this,
m_pfoListName, pPfoList);
86 return STATUS_CODE_SUCCESS;
92 const PfoList *
const pPfoList,
ValidationInfo &validationInfo)
const 94 if (pMCParticleList && pCaloHitList)
120 PfoList allConnectedPfos;
123 PfoList finalStatePfos;
124 for (
const ParticleFlowObject *
const pPfo : allConnectedPfos)
127 finalStatePfos.push_back(pPfo);
149 if (printToScreen && useInterpretedMatching) std::cout <<
"---INTERPRETED-MATCHING-OUTPUT------------------------------------------------------------------" << std::endl;
150 else if (printToScreen) std::cout <<
"---RAW-MATCHING-OUTPUT--------------------------------------------------------------------------" << std::endl;
158 int nNeutrinoPrimaries(0);
159 for (
const MCParticle *
const pMCPrimary : mcPrimaryVector)
162 PfoVector primaryPfoVector;
165 int pfoIndex(0), neutrinoPfoIndex(0);
167 for (
const Pfo *
const pPrimaryPfo : primaryPfoVector)
169 pfoToIdMap.insert(PfoToIdMap::value_type(pPrimaryPfo, ++pfoIndex));
172 if (pRecoNeutrino && !neutrinoPfoToIdMap.count(pRecoNeutrino))
173 neutrinoPfoToIdMap.insert(PfoToIdMap::value_type(pRecoNeutrino, ++neutrinoPfoIndex));
176 PfoSet recoNeutrinos;
177 MCParticleList associatedMCPrimaries;
179 int nCorrectNu(0), nTotalNu(0), nCorrectTB(0), nTotalTB(0), nCorrectCR(0), nTotalCR(0), nFakeNu(0), nFakeCR(0), nSplitNu(0), nSplitCR(0), nLost(0);
180 int mcPrimaryIndex(0), nTargetMatches(0), nTargetNuMatches(0), nTargetCRMatches(0), nTargetGoodNuMatches(0), nTargetNuSplits(0), nTargetNuLosses(0);
181 IntVector mcPrimaryId, mcPrimaryPdg, nMCHitsTotal, nMCHitsU, nMCHitsV, nMCHitsW;
182 FloatVector mcPrimaryE, mcPrimaryPX, mcPrimaryPY, mcPrimaryPZ;
183 FloatVector mcPrimaryVtxX, mcPrimaryVtxY, mcPrimaryVtxZ, mcPrimaryEndX, mcPrimaryEndY, mcPrimaryEndZ;
184 IntVector nPrimaryMatchedPfos, nPrimaryMatchedNuPfos, nPrimaryMatchedCRPfos;
185 IntVector bestMatchPfoId, bestMatchPfoPdg, bestMatchPfoIsRecoNu, bestMatchPfoRecoNuId, bestMatchPfoIsTestBeam;
186 IntVector bestMatchPfoNHitsTotal, bestMatchPfoNHitsU, bestMatchPfoNHitsV, bestMatchPfoNHitsW;
187 IntVector bestMatchPfoNSharedHitsTotal, bestMatchPfoNSharedHitsU, bestMatchPfoNSharedHitsV, bestMatchPfoNSharedHitsW;
189 std::stringstream targetSS;
191 for (
const MCParticle *
const pMCPrimary : mcPrimaryVector)
193 const bool hasMatch(mcToPfoHitSharingMap.count(pMCPrimary) && !mcToPfoHitSharingMap.at(pMCPrimary).empty());
196 if (!hasMatch && !isTargetPrimary)
199 associatedMCPrimaries.push_back(pMCPrimary);
200 const int nTargetPrimaries(associatedMCPrimaries.size());
201 const bool isLastNeutrinoPrimary(++mcPrimaryIndex == nNeutrinoPrimaries);
210 const float targetVertexX(targetVertex.GetX()), targetVertexY(targetVertex.GetY()), targetVertexZ(targetVertex.GetZ());
212 targetSS << (!isTargetPrimary ?
"(Non target) " :
"")
213 <<
"PrimaryId " << mcPrimaryIndex
214 <<
", Nu " << isBeamNeutrinoFinalState
215 <<
", TB " << isBeamParticle
216 <<
", CR " << isCosmicRay
217 <<
", MCPDG " << pMCPrimary->GetParticleId()
218 <<
", Energy " << pMCPrimary->GetEnergy()
219 <<
", Dist. " << (pMCPrimary->GetEndpoint() - pMCPrimary->GetVertex()).GetMagnitude()
220 <<
", nMCHits " << mcPrimaryHitList.size()
225 mcPrimaryId.push_back(mcPrimaryIndex);
226 mcPrimaryPdg.push_back(pMCPrimary->GetParticleId());
227 mcPrimaryE.push_back(pMCPrimary->GetEnergy());
228 mcPrimaryPX.push_back(pMCPrimary->GetMomentum().GetX());
229 mcPrimaryPY.push_back(pMCPrimary->GetMomentum().GetY());
230 mcPrimaryPZ.push_back(pMCPrimary->GetMomentum().GetZ());
231 mcPrimaryVtxX.push_back(pMCPrimary->GetVertex().GetX());
232 mcPrimaryVtxY.push_back(pMCPrimary->GetVertex().GetY());
233 mcPrimaryVtxZ.push_back(pMCPrimary->GetVertex().GetZ());
234 mcPrimaryEndX.push_back(pMCPrimary->GetEndpoint().GetX());
235 mcPrimaryEndY.push_back(pMCPrimary->GetEndpoint().GetY());
236 mcPrimaryEndZ.push_back(pMCPrimary->GetEndpoint().GetZ());
237 nMCHitsTotal.push_back(mcPrimaryHitList.size());
242 int matchIndex(0), nPrimaryMatches(0), nPrimaryNuMatches(0), nPrimaryCRMatches(0), nPrimaryGoodNuMatches(0), nPrimaryNuSplits(0);
248 const CaloHitList &sharedHitList(pfoToSharedHits.second);
249 const CaloHitList &pfoHitList(validationInfo.
GetPfoToHitsMap().at(pfoToSharedHits.first));
253 const bool isGoodMatch(this->
IsGoodMatch(mcPrimaryHitList, pfoHitList, sharedHitList));
255 const int pfoId(pfoToIdMap.at(pfoToSharedHits.first));
258 if (0 == matchIndex++)
260 bestMatchPfoId.push_back(pfoId);
261 bestMatchPfoPdg.push_back(pfoToSharedHits.first->GetParticleId());
262 bestMatchPfoIsRecoNu.push_back(isRecoNeutrinoFinalState ? 1 : 0);
263 bestMatchPfoRecoNuId.push_back(recoNuId);
264 bestMatchPfoIsTestBeam.push_back(isRecoTestBeam ? 1 : 0);
265 bestMatchPfoNHitsTotal.push_back(pfoHitList.size());
269 bestMatchPfoNSharedHitsTotal.push_back(sharedHitList.size());
277 recoVertexX = pRecoVertex->GetPosition().GetX();
278 recoVertexY = pRecoVertex->GetPosition().GetY();
279 recoVertexZ = pRecoVertex->GetPosition().GetZ();
281 catch (
const StatusCodeException &) {}
285 if (isGoodMatch) ++nPrimaryMatches;
287 if (isRecoNeutrinoFinalState)
290 const bool isSplitRecoNeutrino(!recoNeutrinos.empty() && !recoNeutrinos.count(pRecoNeutrino));
291 if (!isSplitRecoNeutrino && isGoodMatch) ++nPrimaryGoodNuMatches;
292 if (isSplitRecoNeutrino && isBeamNeutrinoFinalState && isGoodMatch) ++nPrimaryNuSplits;
293 recoNeutrinos.insert(pRecoNeutrino);
298 if (isRecoNeutrinoFinalState && isGoodMatch) ++nPrimaryNuMatches;
299 if (!isRecoNeutrinoFinalState && isGoodMatch) ++nPrimaryCRMatches;
304 if (isTestBeam && isGoodMatch) ++nPrimaryNuMatches;
305 if (!isTestBeam && isGoodMatch) ++nPrimaryCRMatches;
308 targetSS <<
"-" << (!isGoodMatch ?
"(Below threshold) " :
"")
309 <<
"MatchedPfoId " << pfoId
310 <<
", Nu " << isRecoNeutrinoFinalState;
311 if (isRecoNeutrinoFinalState) targetSS <<
" [NuId: " << recoNuId <<
"]";
312 targetSS <<
", TB " << isRecoTestBeam
313 <<
", CR " << (!isRecoNeutrinoFinalState && !isRecoTestBeam)
314 <<
", PDG " << pfoToSharedHits.first->GetParticleId()
315 <<
", nMatchedHits " << sharedHitList.size()
319 <<
", nPfoHits " << pfoHitList.size()
325 if (mcToPfoHitSharingMap.at(pMCPrimary).empty())
327 targetSS <<
"-No matched Pfo" << std::endl;
328 bestMatchPfoId.push_back(-1); bestMatchPfoPdg.push_back(0); bestMatchPfoIsRecoNu.push_back(0); bestMatchPfoRecoNuId.push_back(-1); bestMatchPfoIsTestBeam.push_back(0);
329 bestMatchPfoNHitsTotal.push_back(0); bestMatchPfoNHitsU.push_back(0); bestMatchPfoNHitsV.push_back(0); bestMatchPfoNHitsW.push_back(0);
330 bestMatchPfoNSharedHitsTotal.push_back(0); bestMatchPfoNSharedHitsU.push_back(0); bestMatchPfoNSharedHitsV.push_back(0); bestMatchPfoNSharedHitsW.push_back(0);
333 nPrimaryMatchedPfos.push_back(nPrimaryMatches);
334 nPrimaryMatchedNuPfos.push_back(nPrimaryNuMatches);
335 nPrimaryMatchedCRPfos.push_back(nPrimaryCRMatches);
336 nTargetMatches += nPrimaryMatches;
337 nTargetNuMatches += nPrimaryNuMatches;
338 nTargetCRMatches += nPrimaryCRMatches;
339 nTargetGoodNuMatches += nPrimaryGoodNuMatches;
340 nTargetNuSplits += nPrimaryNuSplits;
341 if (0 == nPrimaryMatches) ++nTargetNuLosses;
346 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"eventNumber",
m_eventNumber - 1));
347 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcNuanceCode", mcNuanceCode));
348 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isNeutrino", isBeamNeutrinoFinalState));
349 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isBeamParticle", isBeamParticle));
350 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isCosmicRay", isCosmicRay));
351 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"nTargetPrimaries", nTargetPrimaries));
352 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"targetVertexX", targetVertexX));
353 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"targetVertexY", targetVertexY));
354 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"targetVertexZ", targetVertexZ));
355 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"recoVertexX", recoVertexX));
356 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"recoVertexY", recoVertexY));
357 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"recoVertexZ", recoVertexZ));
358 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryId", &mcPrimaryId));
359 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryPdg", &mcPrimaryPdg));
360 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryE", &mcPrimaryE));
361 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryPX", &mcPrimaryPX));
362 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryPY", &mcPrimaryPY));
363 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryPZ", &mcPrimaryPZ));
364 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryVtxX", &mcPrimaryVtxX));
365 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryVtxY", &mcPrimaryVtxY));
366 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryVtxZ", &mcPrimaryVtxZ));
367 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryEndX", &mcPrimaryEndX));
368 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryEndY", &mcPrimaryEndY));
369 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryEndZ", &mcPrimaryEndZ));
370 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryNHitsTotal", &nMCHitsTotal));
371 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryNHitsU", &nMCHitsU));
372 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryNHitsV", &nMCHitsV));
373 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"mcPrimaryNHitsW", &nMCHitsW));
374 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"nPrimaryMatchedPfos", &nPrimaryMatchedPfos));
375 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"nPrimaryMatchedNuPfos", &nPrimaryMatchedNuPfos));
376 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"nPrimaryMatchedCRPfos", &nPrimaryMatchedCRPfos));
377 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoId", &bestMatchPfoId));
378 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoPdg", &bestMatchPfoPdg));
379 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoIsRecoNu", &bestMatchPfoIsRecoNu));
380 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoRecoNuId", &bestMatchPfoRecoNuId));
381 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoNHitsTotal", &bestMatchPfoNHitsTotal));
382 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoNHitsU", &bestMatchPfoNHitsU));
383 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoNHitsV", &bestMatchPfoNHitsV));
384 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoNHitsW", &bestMatchPfoNHitsW));
385 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoNSharedHitsTotal", &bestMatchPfoNSharedHitsTotal));
386 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoNSharedHitsU", &bestMatchPfoNSharedHitsU));
387 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoNSharedHitsV", &bestMatchPfoNSharedHitsV));
388 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoNSharedHitsW", &bestMatchPfoNSharedHitsW));
389 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"nTargetMatches", nTargetMatches));
390 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"nTargetNuMatches", nTargetNuMatches));
391 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"nTargetCRMatches", nTargetCRMatches));
395 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"bestMatchPfoIsTestBeam", &bestMatchPfoIsTestBeam));
400 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"nTargetGoodNuMatches", nTargetGoodNuMatches));
401 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"nTargetNuSplits", nTargetNuSplits));
402 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"nTargetNuLosses", nTargetNuLosses));
406 if (isLastNeutrinoPrimary || isBeamParticle || isCosmicRay)
410 const int interactionTypeInt(static_cast<int>(interactionType));
413 const int isCorrectNu(isBeamNeutrinoFinalState && (nTargetGoodNuMatches == nTargetNuMatches) && (nTargetGoodNuMatches == nTargetPrimaries) && (nTargetCRMatches == 0) && (nTargetNuSplits == 0) && (nTargetNuLosses == 0));
414 const int isCorrectTB(isBeamParticle && (nTargetNuMatches == 1) && (nTargetCRMatches == 0));
415 const int isCorrectCR(isCosmicRay && (nTargetNuMatches == 0) && (nTargetCRMatches == 1));
416 const int isFakeNu(isCosmicRay && (nTargetNuMatches > 0));
417 const int isFakeCR(!isCosmicRay && (nTargetCRMatches > 0));
418 const int isSplitNu(!isCosmicRay && ((nTargetNuMatches > nTargetPrimaries) || (nTargetNuSplits > 0)));
419 const int isSplitCR(isCosmicRay && (nTargetCRMatches > 1));
420 const int isLost(nTargetMatches == 0);
422 std::stringstream outcomeSS;
423 outcomeSS <<
LArInteractionTypeHelper::ToString(interactionType) <<
" (Nuance " << mcNuanceCode <<
", Nu " << isBeamNeutrinoFinalState <<
", TB " << isBeamParticle <<
", CR " << isCosmicRay <<
")" << std::endl;
425 if (isLastNeutrinoPrimary) ++nTotalNu;
426 if (isBeamParticle) ++nTotalTB;
427 if (isCosmicRay) ++nTotalCR;
428 if (isCorrectNu) ++nCorrectNu;
429 if (isCorrectTB) ++nCorrectTB;
430 if (isCorrectCR) ++nCorrectCR;
431 if (isFakeNu) ++nFakeNu;
432 if (isFakeCR) ++nFakeCR;
433 if (isSplitNu) ++nSplitNu;
434 if (isSplitCR) ++nSplitCR;
437 if (isCorrectNu) outcomeSS <<
"IsCorrectNu ";
438 if (isCorrectTB) outcomeSS <<
"IsCorrectTB ";
439 if (isCorrectCR) outcomeSS <<
"IsCorrectCR ";
440 if (isFakeNu) outcomeSS <<
"IsFakeNu ";
441 if (isFakeCR) outcomeSS <<
"IsFakeCR ";
442 if (isSplitNu) outcomeSS <<
"isSplitNu ";
443 if (isSplitCR) outcomeSS <<
"IsSplitCR ";
444 if (isLost) outcomeSS <<
"IsLost ";
445 if (nTargetNuMatches > 0) outcomeSS <<
"(NNuMatches: " << nTargetNuMatches <<
") ";
446 if (nTargetNuLosses > 0) outcomeSS <<
"(NNuLosses: " << nTargetNuLosses <<
") ";
447 if (nTargetNuSplits > 0) outcomeSS <<
"(NNuSplits: " << nTargetNuSplits <<
") ";
448 if (nTargetCRMatches > 0) outcomeSS <<
"(NCRMatches: " << nTargetCRMatches <<
") ";
449 if (printToScreen) std::cout << outcomeSS.str() << std::endl << targetSS.str() << std::endl;
453 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"interactionType", interactionTypeInt));
454 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isCorrectNu", isCorrectNu));
455 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isCorrectTB", isCorrectTB));
456 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isCorrectCR", isCorrectCR));
457 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isFakeNu", isFakeNu));
458 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isFakeCR", isFakeCR));
461 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isSplitNu", isSplitNu));
463 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isSplitCR", isSplitCR));
464 PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(),
m_treeName.c_str(),
"isLost", isLost));
465 PANDORA_MONITORING_API(FillTree(this->GetPandora(),
m_treeName.c_str()));
468 targetSS.str(std::string()); targetSS.clear();
469 recoNeutrinos.clear(); associatedMCPrimaries.clear();
470 nTargetMatches = 0; nTargetNuMatches = 0; nTargetCRMatches = 0; nTargetGoodNuMatches = 0; nTargetNuSplits = 0; nTargetNuLosses = 0;
471 mcPrimaryId.clear(); mcPrimaryPdg.clear(); nMCHitsTotal.clear(); nMCHitsU.clear(); nMCHitsV.clear(); nMCHitsW.clear();
472 mcPrimaryE.clear(); mcPrimaryPX.clear(); mcPrimaryPY.clear(); mcPrimaryPZ.clear();
473 mcPrimaryVtxX.clear(); mcPrimaryVtxY.clear(); mcPrimaryVtxZ.clear(); mcPrimaryEndX.clear(); mcPrimaryEndY.clear(); mcPrimaryEndZ.clear();
474 nPrimaryMatchedPfos.clear(); nPrimaryMatchedNuPfos.clear(); nPrimaryMatchedCRPfos.clear();
475 bestMatchPfoId.clear(); bestMatchPfoPdg.clear(); bestMatchPfoIsRecoNu.clear(); bestMatchPfoRecoNuId.clear(); bestMatchPfoIsTestBeam.clear();
476 bestMatchPfoNHitsTotal.clear(); bestMatchPfoNHitsU.clear(); bestMatchPfoNHitsV.clear(); bestMatchPfoNHitsW.clear();
477 bestMatchPfoNSharedHitsTotal.clear(); bestMatchPfoNSharedHitsU.clear(); bestMatchPfoNSharedHitsV.clear(); bestMatchPfoNSharedHitsW.clear();
481 if (useInterpretedMatching)
483 std::stringstream summarySS;
484 summarySS <<
"---SUMMARY--------------------------------------------------------------------------------------" << std::endl;
485 if (nTotalNu > 0) summarySS <<
"#CorrectNu: " << nCorrectNu <<
"/" << nTotalNu <<
", Fraction: " << (nTotalNu > 0 ?
static_cast<float>(nCorrectNu) / static_cast<float>(nTotalNu) : 0.f) << std::endl;
486 if (nTotalTB > 0) summarySS <<
"#CorrectTB: " << nCorrectTB <<
"/" << nTotalTB <<
", Fraction: " << (nTotalTB > 0 ?
static_cast<float>(nCorrectTB) / static_cast<float>(nTotalTB) : 0.f) << std::endl;
487 if (nTotalCR > 0) summarySS <<
"#CorrectCR: " << nCorrectCR <<
"/" << nTotalCR <<
", Fraction: " << (nTotalCR > 0 ?
static_cast<float>(nCorrectCR) / static_cast<float>(nTotalCR) : 0.f) << std::endl;
488 if (nFakeNu > 0) summarySS <<
"#FakeNu: " << nFakeNu <<
" ";
489 if (nFakeCR > 0) summarySS <<
"#FakeCR: " << nFakeCR <<
" ";
490 if (nSplitNu > 0) summarySS <<
"#SplitNu: " << nSplitNu <<
" ";
491 if (nSplitCR > 0) summarySS <<
"#SplitCR: " << nSplitCR <<
" ";
492 if (nLost > 0) summarySS <<
"#Lost: " << nLost <<
" ";
493 if (nFakeNu || nFakeCR || nSplitNu || nSplitCR || nLost) summarySS << std::endl;
494 if (printToScreen) std::cout << summarySS.str();
497 if (printToScreen) std::cout <<
"------------------------------------------------------------------------------------------------" << std::endl << std::endl;
508 while (this->
GetStrongestPfoMatch(validationInfo, mcPrimaryVector, usedPfos, interpretedMCToPfoHitSharingMap)) {}
512 for (
const MCParticle *
const pMCPrimary : mcPrimaryVector)
516 return ((a.second.size() != b.second.size()) ? a.second.size() > b.second.size() :
LArPfoHelper::SortByNHits(a.first, b.first)); });
525 const MCParticle *pBestMCParticle(
nullptr);
528 for (
const MCParticle *
const pMCPrimary : mcPrimaryVector)
530 if (interpretedMCToPfoHitSharingMap.count(pMCPrimary))
541 if (usedPfos.count(pfoToSharedHits.first))
547 if (pfoToSharedHits.second.size() > bestPfoHitPair.second.size())
549 pBestMCParticle = pMCPrimary;
550 bestPfoHitPair = pfoToSharedHits;
555 if (!pBestMCParticle || !bestPfoHitPair.first)
558 interpretedMCToPfoHitSharingMap[pBestMCParticle].push_back(bestPfoHitPair);
559 usedPfos.insert(bestPfoHitPair.first);
570 for (
const MCParticle *
const pMCPrimary : mcPrimaryVector)
580 if (usedPfos.count(pfoToSharedHits.first))
586 if (pfoToMCParticleHitSharingMap.end() == iter)
588 pfoToMCParticleHitSharingMap[pfoToSharedHits.first].push_back(mcParticleToHits);
592 if (1 != iter->second.size())
593 throw StatusCodeException(STATUS_CODE_FAILURE);
597 if (mcParticleToHits.second.size() > originalMCParticleToHits.second.size())
598 originalMCParticleToHits = mcParticleToHits;
603 for (
const auto &mapEntry : pfoToMCParticleHitSharingMap)
614 const float purity((recoHits.size() > 0) ? static_cast<float>(sharedHits.size()) / static_cast<float>(recoHits.size()) : 0.
f);
615 const float completeness((trueHits.size() > 0) ? static_cast<float>(sharedHits.size()) / static_cast<float>(trueHits.size()) : 0.
f);
624 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"CaloHitListName",
m_caloHitListName));
625 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"MCParticleListName",
m_mcParticleListName));
626 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"PfoListName",
m_pfoListName));
628 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
631 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
634 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
637 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
640 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
643 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
646 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
649 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
652 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
655 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
658 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
661 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
666 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"OutputTree",
m_treeName));
667 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"OutputFile",
m_fileName));
669 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
673 return STATUS_CODE_SUCCESS;
static bool SortByNHits(const pandora::ParticleFlowObject *const pLhs, const pandora::ParticleFlowObject *const pRhs)
Sort pfos by number of constituent hits.
Header file for the pfo helper class.
bool m_selectInputHits
whether to select input hits
~EventValidationAlgorithm()
Destructor.
Header file for the interaction type helper class.
pandora::StatusCode Run()
std::string m_mcParticleListName
Name of input MC particle list.
std::string m_pfoListName
Name of input Pfo list.
void SetAllMCParticleToHitsMap(const LArMCParticleHelper::MCContributionMap &allMCParticleToHitsMap)
Set the all mc particle to hits map.
bool m_useTrueNeutrinosOnly
Whether to consider only mc particles that were neutrino induced.
int m_eventNumber
The event number.
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.
bool GetStrongestPfoMatch(const ValidationInfo &validationInfo, const pandora::MCParticleVector &mcPrimaryVector, pandora::PfoSet &usedPfos, LArMCParticleHelper::MCParticleToPfoHitSharingMap &interpretedMCToPfoHitSharingMap) const
Get the strongest pfo match (most matched hits) between an available mc primary and an available pfo...
static void GetPfoMCParticleHitSharingMaps(const PfoContributionMap &pfoToReconstructable2DHitsMap, const MCContributionMapVector &selectedMCParticleToHitsMaps, PfoToMCParticleHitSharingMap &pfoToMCParticleHitSharingMap, MCParticleToPfoHitSharingMap &mcParticleToPfoHitSharingMap)
Get the mappings from Pfo -> pair (reconstructable MCparticles, number of reconstructable 2D hits sha...
void SetTargetMCParticleToHitsMap(const LArMCParticleHelper::MCContributionMap &targetMCParticleToHitsMap)
Set the target mc particle to hits map.
std::vector< PfoCaloHitListPair > PfoToSharedHitsVector
void GetRemainingPfoMatches(const ValidationInfo &validationInfo, const pandora::MCParticleVector &mcPrimaryVector, const pandora::PfoSet &usedPfos, LArMCParticleHelper::MCParticleToPfoHitSharingMap &interpretedMCToPfoHitSharingMap) const
Get the best matches for any pfos left-over after the strong matching procedure.
bool m_useSmallPrimaries
Whether to consider matches to mc primaries with fewer than m_matchingMinPrimaryHits.
std::vector< int > IntVector
bool m_printAllToScreen
Whether to print all/raw matching details to screen.
float m_matchingMinCompleteness
The minimum particle completeness to declare a match.
Header file for the lar monitoring helper helper class.
bool m_testBeamMode
Whether pandora is reconstructing test beam particles.
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
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
float m_maxPhotonPropagation
the maximum photon propagation length
std::unordered_map< const pandora::ParticleFlowObject *, unsigned int > PfoToIdMap
bool m_printMatchingToScreen
Whether to print matching output to screen.
static void GetOrderedMCParticleVector(const LArMCParticleHelper::MCContributionMapVector &selectedMCParticleToGoodHitsMaps, pandora::MCParticleVector &orderedMCParticleVector)
Order input MCParticles by their number of hits.
void PrintAllMatches(const ValidationInfo &validationInfo) const
Print all/raw matching information to screen.
static bool IsCosmicRay(const pandora::MCParticle *const pMCParticle)
Return true if passed a primary cosmic ray MCParticle.
Header file for the event validation algorithm.
static void SelectReconstructableMCParticles(const pandora::MCParticleList *pMCParticleList, const pandora::CaloHitList *pCaloHitList, const PrimaryParameters ¶meters, std::function< bool(const pandora::MCParticle *const)> fCriteria, MCContributionMap &selectedMCParticlesToHitsMap)
Select primary, reconstructable mc particles that match given criteria.
float m_minHitSharingFraction
Minimum fraction of energy deposited by selected primary in a single "good" hit.
bool m_writeToTree
Whether to write all/raw matching details to tree.
std::string m_treeName
Name of output tree.
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...
static bool IsTestBeam(const pandora::ParticleFlowObject *const pPfo)
Whether a pfo is a test beam particle.
std::unordered_map< const pandora::MCParticle *, pandora::CaloHitList > MCContributionMap
static void GetPfoToReconstructable2DHitsMap(const pandora::PfoList &pfoList, const MCContributionMap &selectedMCParticleToHitsMap, PfoContributionMap &pfoToReconstructable2DHitsMap)
Get mapping from Pfo to reconstructable 2D hits (=good hits belonging to a selected reconstructable M...
void SetPfoToHitsMap(const LArMCParticleHelper::PfoContributionMap &pfoToHitsMap)
Set the pfo to hits map.
static std::string ToString(const InteractionType interactionType)
Get a string representation of an interaction type.
static bool IsBeamParticle(const pandora::MCParticle *const pMCParticle)
Returns true if passed a primary beam MCParticle.
const LArMCParticleHelper::PfoContributionMap & GetPfoToHitsMap() const
Get the pfo to hits map.
std::map< const pandora::MCParticle *, PfoToSharedHitsVector > MCParticleToPfoHitSharingMap
void PrintInterpretedMatches(const ValidationInfo &validationInfo) const
Print interpreted matching information to screen.
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.
float m_minHitSharingFraction
the minimum Hit sharing fraction
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.
bool m_selectInputHits
Whether to use only hits passing mc-based quality (is "reconstructable") checks.
void SetMCToPfoHitSharingMap(const LArMCParticleHelper::MCParticleToPfoHitSharingMap &mcToPfoHitSharingMap)
Set the mc to pfo hit sharing map.
int m_fileIdentifier
The input file identifier.
InteractionType
InteractionType enum.
void SetInterpretedMCToPfoHitSharingMap(const LArMCParticleHelper::MCParticleToPfoHitSharingMap &interpretedMCToPfoHitSharingMap)
Set the interpreted mc to pfo hit sharing map.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
std::map< const pandora::ParticleFlowObject *, MCParticleToSharedHitsVector > PfoToMCParticleHitSharingMap
std::string m_fileName
Name of output file.
const LArMCParticleHelper::MCContributionMap & GetAllMCParticleToHitsMap() const
Get the all mc particle to hits map.
std::pair< const pandora::MCParticle *, pandora::CaloHitList > MCParticleCaloHitListPair
static bool IsFinalState(const pandora::ParticleFlowObject *const pPfo)
Whether a pfo is a primary parent particle.
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.
float m_matchingMinPurity
The minimum particle purity to declare a match.
unsigned int m_matchingMinSharedHits
The minimum number of shared hits used in matching scheme.
static bool IsNeutrinoFinalState(const pandora::ParticleFlowObject *const pPfo)
Whether a pfo is a final-state particle from a neutrino (or antineutrino) interaction.
float m_maxPhotonPropagation
Maximum distance travelled by photon, downstream of a track, in mc particle hierarchy.
std::unordered_map< const pandora::ParticleFlowObject *, pandora::CaloHitList > PfoContributionMap
void WriteInterpretedMatches(const ValidationInfo &validationInfo) const
Write interpreted matching information to tree.
const LArMCParticleHelper::MCParticleToPfoHitSharingMap & GetInterpretedMCToPfoHitSharingMap() const
Get the interpreted mc to pfo hit sharing map.
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...
std::string m_caloHitListName
Name of input calo hit list.
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.
const LArMCParticleHelper::MCParticleToPfoHitSharingMap & GetMCToPfoHitSharingMap() const
Get the mc to pfo hit sharing map.
void InterpretMatching(const ValidationInfo &validationInfo, LArMCParticleHelper::MCParticleToPfoHitSharingMap &interpretedMCToPfoHitSharingMap) const
Apply an interpretative matching procedure to the comprehensive matches in the provided validation in...
const LArMCParticleHelper::MCContributionMap & GetTargetMCParticleToHitsMap() const
Get the target mc particle to hits map.
std::pair< const pandora::ParticleFlowObject *, pandora::CaloHitList > PfoCaloHitListPair