20 TChain *pTChain =
new TChain(
"Validation",
"pTChain");
21 pTChain->Add(inputFiles.c_str());
26 int nEvents(0), nProcessedEvents(0);
27 const int nChainEntries(pTChain->GetEntries());
29 for (
int iEntry = 0; iEntry < nChainEntries; )
32 iEntry +=
ReadNextEvent(pTChain, iEntry, simpleMCEvent, parameters);
37 if (nEvents % 50 == 0)
38 std::cout <<
"nEvents " << nEvents <<
"\r" << std::flush;
46 CountPfoMatches(simpleMCEvent, parameters, interactionCountingMap, interactionTargetResultMap);
57 int thisEventNumber(0), iTarget(0);
58 const int nChainEntries(pTChain->GetEntries());
60 pTChain->SetBranchAddress(
"eventNumber", &thisEventNumber);
61 pTChain->SetBranchAddress(
"fileIdentifier", &simpleMCEvent.
m_fileIdentifier);
62 pTChain->GetEntry(iEntry);
65 while (iEntry + iTarget < nChainEntries)
70 pTChain->SetBranchAddress(
"mcNuanceCode", &simpleMCTarget.
m_mcNuanceCode);
71 pTChain->SetBranchAddress(
"isCosmicRay", &simpleMCTarget.
m_isCosmicRay);
75 pTChain->SetBranchAddress(
"recoVertexX", &simpleMCTarget.
m_recoVertex.
m_x);
76 pTChain->SetBranchAddress(
"recoVertexY", &simpleMCTarget.
m_recoVertex.
m_y);
77 pTChain->SetBranchAddress(
"recoVertexZ", &simpleMCTarget.
m_recoVertex.
m_z);
78 pTChain->SetBranchAddress(
"isCorrectCR", &simpleMCTarget.
m_isCorrectCR);
79 pTChain->SetBranchAddress(
"isFakeCR", &simpleMCTarget.
m_isFakeCR);
80 pTChain->SetBranchAddress(
"isSplitCR", &simpleMCTarget.
m_isSplitCR);
81 pTChain->SetBranchAddress(
"isLost", &simpleMCTarget.
m_isLost);
82 pTChain->SetBranchAddress(
"nTargetMatches", &simpleMCTarget.
m_nTargetMatches);
88 pTChain->SetBranchAddress(
"isBeamParticle", &simpleMCTarget.
m_isBeamParticle);
89 pTChain->SetBranchAddress(
"isCorrectTB", &simpleMCTarget.
m_isCorrectTB);
93 pTChain->SetBranchAddress(
"isNeutrino", &simpleMCTarget.
m_isNeutrino);
94 pTChain->SetBranchAddress(
"isCorrectNu", &simpleMCTarget.
m_isCorrectNu);
95 pTChain->SetBranchAddress(
"isFakeNu", &simpleMCTarget.
m_isFakeNu);
96 pTChain->SetBranchAddress(
"isSplitNu", &simpleMCTarget.
m_isSplitNu);
100 IntVector *pMCPrimaryId(
nullptr), *pMCPrimaryPdg(
nullptr), *pNMCHitsTotal(
nullptr), *pNMCHitsU(
nullptr), *pNMCHitsV(
nullptr), *pNMCHitsW(
nullptr);
101 FloatVector *pMCPrimaryE(
nullptr), *pMCPrimaryPX(
nullptr), *pMCPrimaryPY(
nullptr), *pMCPrimaryPZ(
nullptr);
102 FloatVector *pMCPrimaryVtxX(
nullptr), *pMCPrimaryVtxY(
nullptr), *pMCPrimaryVtxZ(
nullptr), *pMCPrimaryEndX(
nullptr), *pMCPrimaryEndY(
nullptr), *pMCPrimaryEndZ(
nullptr);
103 IntVector *pNPrimaryMatchedPfos(
nullptr), *pNPrimaryMatchedNuPfos(
nullptr), *pNPrimaryMatchedCRPfos(
nullptr);
104 IntVector *pBestMatchPfoId(
nullptr), *pBestMatchPfoPdg(
nullptr), *pBestMatchPfoIsRecoNu(
nullptr), *pBestMatchPfoRecoNuId(
nullptr), *pBestMatchPfoIsTestBeam(
nullptr);
105 IntVector *pBestMatchPfoNHitsTotal(
nullptr), *pBestMatchPfoNHitsU(
nullptr), *pBestMatchPfoNHitsV(
nullptr), *pBestMatchPfoNHitsW(
nullptr);
106 IntVector *pBestMatchPfoNSharedHitsTotal(
nullptr), *pBestMatchPfoNSharedHitsU(
nullptr), *pBestMatchPfoNSharedHitsV(
nullptr), *pBestMatchPfoNSharedHitsW(
nullptr);
108 pTChain->SetBranchAddress(
"mcPrimaryId", &pMCPrimaryId);
109 pTChain->SetBranchAddress(
"mcPrimaryPdg", &pMCPrimaryPdg);
110 pTChain->SetBranchAddress(
"mcPrimaryE", &pMCPrimaryE);
111 pTChain->SetBranchAddress(
"mcPrimaryPX", &pMCPrimaryPX);
112 pTChain->SetBranchAddress(
"mcPrimaryPY", &pMCPrimaryPY);
113 pTChain->SetBranchAddress(
"mcPrimaryPZ", &pMCPrimaryPZ);
114 pTChain->SetBranchAddress(
"mcPrimaryVtxX", &pMCPrimaryVtxX);
115 pTChain->SetBranchAddress(
"mcPrimaryVtxY", &pMCPrimaryVtxY);
116 pTChain->SetBranchAddress(
"mcPrimaryVtxZ", &pMCPrimaryVtxZ);
117 pTChain->SetBranchAddress(
"mcPrimaryEndX", &pMCPrimaryEndX);
118 pTChain->SetBranchAddress(
"mcPrimaryEndY", &pMCPrimaryEndY);
119 pTChain->SetBranchAddress(
"mcPrimaryEndZ", &pMCPrimaryEndZ);
120 pTChain->SetBranchAddress(
"mcPrimaryNHitsTotal", &pNMCHitsTotal);
121 pTChain->SetBranchAddress(
"mcPrimaryNHitsU", &pNMCHitsU);
122 pTChain->SetBranchAddress(
"mcPrimaryNHitsV", &pNMCHitsV);
123 pTChain->SetBranchAddress(
"mcPrimaryNHitsW", &pNMCHitsW);
124 pTChain->SetBranchAddress(
"nPrimaryMatchedPfos", &pNPrimaryMatchedPfos);
125 pTChain->SetBranchAddress(
"nPrimaryMatchedCRPfos", &pNPrimaryMatchedCRPfos);
126 pTChain->SetBranchAddress(
"bestMatchPfoNHitsTotal", &pBestMatchPfoNHitsTotal);
127 pTChain->SetBranchAddress(
"bestMatchPfoNHitsU", &pBestMatchPfoNHitsU);
128 pTChain->SetBranchAddress(
"bestMatchPfoNHitsV", &pBestMatchPfoNHitsV);
129 pTChain->SetBranchAddress(
"bestMatchPfoNHitsW", &pBestMatchPfoNHitsW);
130 pTChain->SetBranchAddress(
"bestMatchPfoId", &pBestMatchPfoId);
131 pTChain->SetBranchAddress(
"bestMatchPfoPdg", &pBestMatchPfoPdg);
132 pTChain->SetBranchAddress(
"bestMatchPfoNSharedHitsTotal", &pBestMatchPfoNSharedHitsTotal);
133 pTChain->SetBranchAddress(
"bestMatchPfoNSharedHitsU", &pBestMatchPfoNSharedHitsU);
134 pTChain->SetBranchAddress(
"bestMatchPfoNSharedHitsV", &pBestMatchPfoNSharedHitsV);
135 pTChain->SetBranchAddress(
"bestMatchPfoNSharedHitsW", &pBestMatchPfoNSharedHitsW);
139 pTChain->SetBranchAddress(
"bestMatchPfoIsTB", &pBestMatchPfoIsTestBeam);
143 pTChain->SetBranchAddress(
"nPrimaryMatchedNuPfos", &pNPrimaryMatchedNuPfos);
144 pTChain->SetBranchAddress(
"bestMatchPfoIsRecoNu", &pBestMatchPfoIsRecoNu);
145 pTChain->SetBranchAddress(
"bestMatchPfoRecoNuId", &pBestMatchPfoRecoNuId);
151 pTChain->GetEntry(iEntry + iTarget++);
159 simpleMCPrimary.
m_primaryId = pMCPrimaryId->at(iPrimary);
160 simpleMCPrimary.
m_pdgCode = pMCPrimaryPdg->at(iPrimary);
161 simpleMCPrimary.
m_energy = pMCPrimaryE->at(iPrimary);
165 simpleMCPrimary.
m_vertex.
m_x = pMCPrimaryVtxX->at(iPrimary);
166 simpleMCPrimary.
m_vertex.
m_y = pMCPrimaryVtxY->at(iPrimary);
167 simpleMCPrimary.
m_vertex.
m_z = pMCPrimaryVtxZ->at(iPrimary);
168 simpleMCPrimary.
m_endpoint.
m_x = pMCPrimaryEndX->at(iPrimary);
169 simpleMCPrimary.
m_endpoint.
m_y = pMCPrimaryEndY->at(iPrimary);
170 simpleMCPrimary.
m_endpoint.
m_z = pMCPrimaryEndZ->at(iPrimary);
172 simpleMCPrimary.
m_nMCHitsU = pNMCHitsU->at(iPrimary);
173 simpleMCPrimary.
m_nMCHitsV = pNMCHitsV->at(iPrimary);
174 simpleMCPrimary.
m_nMCHitsW = pNMCHitsW->at(iPrimary);
206 pTChain->ResetBranchAddresses();
214 std::cout <<
"---INTERPRETED-MATCHING-OUTPUT------------------------------------------------------------------" << std::endl;
217 int nCorrectNu(0), nTotalNu(0), nCorrectTB(0), nTotalTB(0), nCorrectCR(0), nTotalCR(0), nFakeNu(0), nFakeCR(0), nSplitNu(0), nSplitCR(0), nLost(0);
221 std::cout << std::endl << ToString(static_cast<InteractionType>(simpleMCTarget.m_interactionType))
222 <<
" (Nuance " << simpleMCTarget.m_mcNuanceCode <<
", Nu " << simpleMCTarget.m_isNeutrino;
223 if (!
PassFiducialCut(simpleMCTarget, parameters) && simpleMCTarget.m_isNeutrino) std::cout <<
" [NonFid]";
224 std::cout <<
", TB " << simpleMCTarget.m_isBeamParticle <<
", CR " << simpleMCTarget.m_isCosmicRay <<
")" << std::endl;
226 std::stringstream ss;
227 if (simpleMCTarget.m_isCorrectNu) ss <<
"IsCorrectNu ";
228 if (simpleMCTarget.m_isCorrectTB) ss <<
"IsCorrectTB ";
229 if (simpleMCTarget.m_isCorrectCR) ss <<
"IsCorrectCR ";
230 if (simpleMCTarget.m_isFakeNu) ss <<
"IsFakeNu ";
231 if (simpleMCTarget.m_isFakeCR) ss <<
"IsFakeCR ";
232 if (simpleMCTarget.m_isSplitNu) ss <<
"IsSplitNu ";
233 if (simpleMCTarget.m_isSplitCR) ss <<
"IsSplitCR ";
234 if (simpleMCTarget.m_isLost) ss <<
"IsLost ";
235 if (simpleMCTarget.m_nTargetNuMatches > 0) ss <<
"(NNuMatches: " << simpleMCTarget.m_nTargetNuMatches <<
") ";
236 if (simpleMCTarget.m_nTargetNuSplits > 0) ss <<
"(NNuSplits: " << simpleMCTarget.m_nTargetNuSplits <<
") ";
237 if (simpleMCTarget.m_nTargetNuLosses > 0) ss <<
"(NNuLosses: " << simpleMCTarget.m_nTargetNuLosses <<
") ";
238 if (simpleMCTarget.m_nTargetCRMatches > 0) ss <<
"(NCRMatches: " << simpleMCTarget.m_nTargetCRMatches <<
") ";
239 std::cout << ss.str() << std::endl;
241 if (simpleMCTarget.m_isNeutrino) ++nTotalNu;
242 if (simpleMCTarget.m_isBeamParticle) ++nTotalTB;
243 if (simpleMCTarget.m_isCosmicRay) ++nTotalCR;
244 if (simpleMCTarget.m_isCorrectNu) ++nCorrectNu;
245 if (simpleMCTarget.m_isCorrectTB) ++nCorrectTB;
246 if (simpleMCTarget.m_isCorrectCR) ++nCorrectCR;
247 if (simpleMCTarget.m_isFakeNu) ++nFakeNu;
248 if (simpleMCTarget.m_isFakeCR) ++nFakeCR;
249 if (simpleMCTarget.m_isSplitNu) ++nSplitNu;
250 if (simpleMCTarget.m_isSplitCR) ++nSplitCR;
251 if (simpleMCTarget.m_isLost) ++nLost;
253 for (
const SimpleMCPrimary &simpleMCPrimary : simpleMCTarget.m_mcPrimaryList)
255 std::cout <<
"PrimaryId " << simpleMCPrimary.m_primaryId
256 <<
", Nu " << simpleMCTarget.m_isNeutrino
257 <<
", TB " << simpleMCTarget.m_isBeamParticle
258 <<
", CR " << simpleMCTarget.m_isCosmicRay
259 <<
", MCPDG " << simpleMCPrimary.m_pdgCode
260 <<
", Energy " << simpleMCPrimary.m_energy
261 <<
", Dist. " << std::sqrt(
262 (simpleMCPrimary.m_vertex.m_x - simpleMCPrimary.m_endpoint.m_x) * (simpleMCPrimary.m_vertex.m_x - simpleMCPrimary.m_endpoint.m_x) +
263 (simpleMCPrimary.m_vertex.m_y - simpleMCPrimary.m_endpoint.m_y) * (simpleMCPrimary.m_vertex.m_y - simpleMCPrimary.m_endpoint.m_y) +
264 (simpleMCPrimary.m_vertex.m_z - simpleMCPrimary.m_endpoint.m_z) * (simpleMCPrimary.m_vertex.m_z - simpleMCPrimary.m_endpoint.m_z))
265 <<
", nMCHits " << simpleMCPrimary.m_nMCHitsTotal
266 <<
" (" << simpleMCPrimary.m_nMCHitsU
267 <<
", " << simpleMCPrimary.m_nMCHitsV
268 <<
", " << simpleMCPrimary.m_nMCHitsW <<
")" << std::endl;
270 if (0 == simpleMCPrimary.m_nPrimaryMatchedPfos)
272 std::cout <<
"-No matched Pfo" << std::endl;
276 std::cout <<
"-MatchedPfoId " << simpleMCPrimary.m_bestMatchPfoId;
277 if (simpleMCPrimary.m_nPrimaryMatchedPfos > 1) std::cout <<
" (NMatches " << simpleMCPrimary.m_nPrimaryMatchedPfos <<
")";
278 std::cout <<
", Nu " << simpleMCPrimary.m_bestMatchPfoIsRecoNu;
279 if (simpleMCPrimary.m_bestMatchPfoIsRecoNu) std::cout <<
" [NuId: " << simpleMCPrimary.m_bestMatchPfoRecoNuId <<
"]";
280 std::cout <<
", TB " << (simpleMCPrimary.m_bestMatchPfoIsTestBeam)
281 <<
", CR " << (!simpleMCPrimary.m_bestMatchPfoIsRecoNu && !simpleMCPrimary.m_bestMatchPfoIsTestBeam)
282 <<
", PDG " << simpleMCPrimary.m_bestMatchPfoPdgCode
283 <<
", nMatchedHits " << simpleMCPrimary.m_bestMatchPfoNSharedHitsTotal
284 <<
" (" << simpleMCPrimary.m_bestMatchPfoNSharedHitsU
285 <<
", " << simpleMCPrimary.m_bestMatchPfoNSharedHitsV
286 <<
", " << simpleMCPrimary.m_bestMatchPfoNSharedHitsW <<
")"
287 <<
", nPfoHits " << simpleMCPrimary.m_bestMatchPfoNHitsTotal
288 <<
" (" << simpleMCPrimary.m_bestMatchPfoNHitsU
289 <<
", " << simpleMCPrimary.m_bestMatchPfoNHitsV
290 <<
", " << simpleMCPrimary.m_bestMatchPfoNHitsW <<
")" << std::endl;
294 std::stringstream summarySS;
295 summarySS << std::endl <<
"---SUMMARY--------------------------------------------------------------------------------------" << std::endl;
296 if (nTotalNu > 0) summarySS <<
"#CorrectNu: " << nCorrectNu <<
"/" << nTotalNu <<
", Fraction: " << (nTotalNu > 0 ?
static_cast<float>(nCorrectNu) /
static_cast<float>(nTotalNu) : 0.f) << std::endl;
297 if (nTotalTB > 0) summarySS <<
"#CorrectTB: " << nCorrectTB <<
"/" << nTotalTB <<
", Fraction: " << (nTotalTB > 0 ?
static_cast<float>(nCorrectTB) /
static_cast<float>(nTotalTB) : 0.f) << std::endl;
298 if (nTotalCR > 0) summarySS <<
"#CorrectCR: " << nCorrectCR <<
"/" << nTotalCR <<
", Fraction: " << (nTotalCR > 0 ?
static_cast<float>(nCorrectCR) /
static_cast<float>(nTotalCR) : 0.f) << std::endl;
299 if (nFakeNu > 0) summarySS <<
"#FakeNu: " << nFakeNu <<
" ";
300 if (nFakeCR > 0) summarySS <<
"#FakeCR: " << nFakeCR <<
" ";
301 if (nSplitNu > 0) summarySS <<
"#SplitNu: " << nSplitNu <<
" ";
302 if (nSplitCR > 0) summarySS <<
"#SplitCR: " << nSplitCR <<
" ";
303 if (nLost > 0) summarySS <<
"#Lost: " << nLost <<
" ";
304 if (nFakeNu || nFakeCR || nSplitNu || nSplitCR || nLost) summarySS << std::endl;
305 std::cout << summarySS.str();
306 std::cout <<
"------------------------------------------------------------------------------------------------" << std::endl;
316 if ((!
PassFiducialCut(simpleMCTarget, parameters) && simpleMCTarget.m_isNeutrino) ||
317 (parameters.
m_triggeredBeamOnly && simpleMCTarget.m_isBeamParticle && simpleMCTarget.m_mcNuanceCode != 2001))
323 targetResult.
m_isCorrect = (simpleMCTarget.m_isNeutrino && simpleMCTarget.m_isCorrectNu) ||
324 (simpleMCTarget.m_isBeamParticle && simpleMCTarget.m_isCorrectTB) ||
325 (simpleMCTarget.m_isCosmicRay && simpleMCTarget.m_isCorrectCR);
327 if (simpleMCTarget.m_nTargetMatches > 0)
330 targetResult.
m_vertexOffset = simpleMCTarget.m_recoVertex - simpleMCTarget.m_targetVertex;
336 for (
const SimpleMCPrimary &simpleMCPrimary : simpleMCTarget.m_mcPrimaryList)
341 CountingDetails &countingDetails = interactionCountingMap[interactionType][expectedPrimary];
345 bool incorrectMatchToCR(parameters.
m_testBeamMode ? (simpleMCTarget.m_isCosmicRay == simpleMCPrimary.m_bestMatchPfoIsTestBeam) : (simpleMCTarget.m_isCosmicRay == simpleMCPrimary.m_bestMatchPfoIsRecoNu));
347 if ((simpleMCPrimary.m_bestMatchPfoId >= 0) && incorrectMatchToCR)
353 if (0 == simpleMCPrimary.m_nPrimaryMatchedPfos) ++countingDetails.
m_nMatch0;
354 else if (1 == simpleMCPrimary.m_nPrimaryMatchedPfos) ++countingDetails.
m_nMatch1;
355 else if (2 == simpleMCPrimary.m_nPrimaryMatchedPfos) ++countingDetails.
m_nMatch2;
358 primaryResult.
m_nPfoMatches = simpleMCPrimary.m_nPrimaryMatchedPfos;
362 primaryResult.
m_bestMatchCompleteness = (simpleMCPrimary.m_nMCHitsTotal > 0) ?
static_cast<float>(simpleMCPrimary.m_bestMatchPfoNSharedHitsTotal) /
static_cast<float>(simpleMCPrimary.m_nMCHitsTotal) : 0.f;
363 primaryResult.
m_bestMatchPurity = (simpleMCPrimary.m_bestMatchPfoNHitsTotal > 0) ?
static_cast<float>(simpleMCPrimary.m_bestMatchPfoNSharedHitsTotal) /
static_cast<float>(simpleMCPrimary.m_bestMatchPfoNHitsTotal) : 0.f;
372 const float pTot(std::sqrt(simpleMCPrimary.m_momentum.m_x * simpleMCPrimary.m_momentum.m_x + simpleMCPrimary.m_momentum.m_y * simpleMCPrimary.m_momentum.m_y + simpleMCPrimary.m_momentum.m_z * simpleMCPrimary.m_momentum.m_z));
376 interactionTargetResultMap[interactionType].push_back(targetResult);
385 throw std::invalid_argument(
"Parameters has fiducial cuts for uBooNE and SBND");
400 const float eVx(256.35), eVy(233.), eVz(1036.8);
401 const float xBorder(10.), yBorder(20.), zBorder(10.);
417 const float eVx(400.f), eVy(400.f), eVz(500.f);
418 const float xBorder(10.f), yBorder(20.f), zBorder(10.f);
436 unsigned int nMuons(0), nElectrons(0), nProtons(0), nPiPlus(0), nPiMinus(0), nNeutrons(0), nPhotons(0);
438 for (
const SimpleMCPrimary &simpleMCPrimaryInList : simpleMCPrimaryList)
440 if (&simpleMCPrimary == &simpleMCPrimaryInList)
442 if ((0 == nMuons) && (13 == std::fabs(simpleMCPrimaryInList.m_pdgCode)))
return MUON;
443 if ((0 == nElectrons) && (11 == std::fabs(simpleMCPrimaryInList.m_pdgCode)))
return ELECTRON;
444 if ((0 == nProtons) && (2212 == std::fabs(simpleMCPrimaryInList.m_pdgCode)))
return PROTON1;
445 if ((1 == nProtons) && (2212 == std::fabs(simpleMCPrimaryInList.m_pdgCode)))
return PROTON2;
446 if ((2 == nProtons) && (2212 == std::fabs(simpleMCPrimaryInList.m_pdgCode)))
return PROTON3;
447 if ((3 == nProtons) && (2212 == std::fabs(simpleMCPrimaryInList.m_pdgCode)))
return PROTON4;
448 if ((4 == nProtons) && (2212 == std::fabs(simpleMCPrimaryInList.m_pdgCode)))
return PROTON5;
449 if ((0 == nPiPlus) && (211 == simpleMCPrimaryInList.m_pdgCode))
return PIPLUS;
450 if ((0 == nPiMinus) && (-211 == simpleMCPrimaryInList.m_pdgCode))
return PIMINUS;
451 if ((0 == nPhotons) && (22 == simpleMCPrimaryInList.m_pdgCode))
return PHOTON1;
452 if ((1 == nPhotons) && (22 == simpleMCPrimaryInList.m_pdgCode))
return PHOTON2;
455 if (13 == std::fabs(simpleMCPrimaryInList.m_pdgCode)) ++nMuons;
456 else if (11 == std::fabs(simpleMCPrimaryInList.m_pdgCode)) ++nElectrons;
457 else if (2212 == std::fabs(simpleMCPrimaryInList.m_pdgCode)) ++nProtons;
458 else if (211 == simpleMCPrimaryInList.m_pdgCode) ++nPiPlus;
459 else if (-211 == simpleMCPrimaryInList.m_pdgCode) ++nPiMinus;
460 else if (2112 == std::fabs(simpleMCPrimaryInList.m_pdgCode)) ++nNeutrons;
461 else if (22 == simpleMCPrimaryInList.m_pdgCode) ++nPhotons;
471 const unsigned int absMCPdgCode(std::fabs(simpleMCPrimary.
m_pdgCode));
473 if (((absMCPdgCode == 13 || absMCPdgCode == 2212 || absMCPdgCode == 211) && (13 != std::fabs(bestMatchPfoPdgCode) && 211 != std::fabs(bestMatchPfoPdgCode))) ||
474 ((absMCPdgCode == 22 || absMCPdgCode == 11) && (11 != std::fabs(bestMatchPfoPdgCode))) )
486 std::cout << std::fixed;
487 std::cout << std::setprecision(1);
489 std::ofstream mapFile;
492 for (
const InteractionCountingMap::value_type &interactionTypeMapEntry : interactionCountingMap)
495 const CountingMap &countingMap(interactionTypeMapEntry.second);
496 std::cout << std::endl <<
ToString(interactionType) << std::endl;
499 mapFile << std::endl <<
ToString(interactionType) << std::endl;
501 for (
const CountingMap::value_type &countingMapEntry : countingMap)
506 std::cout <<
"-" <<
ToString(expectedPrimary) <<
": nEvents: " << countingDetails.
m_nTotal
507 <<
", nPfos |0: " << ((countingDetails.
m_nTotal > 0) ? 100.f *
static_cast<float>(countingDetails.
m_nMatch0) /
static_cast<float>(countingDetails.
m_nTotal) : 0.f)
508 <<
"%|, |1: " << ((countingDetails.
m_nTotal > 0) ? 100.f *
static_cast<float>(countingDetails.
m_nMatch1) /
static_cast<float>(countingDetails.
m_nTotal) : 0.f)
509 <<
"%|, |2: " << ((countingDetails.
m_nTotal > 0) ? 100.f *
static_cast<float>(countingDetails.
m_nMatch2) /
static_cast<float>(countingDetails.
m_nTotal) : 0.f)
510 <<
"%|, |3+: " << ((countingDetails.
m_nTotal > 0) ? 100.f *
static_cast<float>(countingDetails.
m_nMatch3Plus) /
static_cast<float>(countingDetails.
m_nTotal) : 0.f)
511 <<
"%|, correctId " << ((countingDetails.
m_nTotal - countingDetails.
m_nMatch0 > 0) ? 100.f *
static_cast<float>(countingDetails.
m_correctId) /
static_cast<float>(countingDetails.
m_nTotal - countingDetails.
m_nMatch0) : 0.f)
516 mapFile <<
"-" <<
ToString(expectedPrimary) <<
": nEvents: " << countingDetails.
m_nTotal
517 <<
", nPfos |0: " << ((countingDetails.
m_nTotal > 0) ? 100.f *
static_cast<float>(countingDetails.
m_nMatch0) /
static_cast<float>(countingDetails.
m_nTotal) : 0.f)
518 <<
"%|, |1: " << ((countingDetails.
m_nTotal > 0) ? 100.f *
static_cast<float>(countingDetails.
m_nMatch1) /
static_cast<float>(countingDetails.
m_nTotal) : 0.f)
519 <<
"%|, |2: " << ((countingDetails.
m_nTotal > 0) ? 100.f *
static_cast<float>(countingDetails.
m_nMatch2) /
static_cast<float>(countingDetails.
m_nTotal) : 0.f)
520 <<
"%|, |3+: " << ((countingDetails.
m_nTotal > 0) ? 100.f *
static_cast<float>(countingDetails.
m_nMatch3Plus) /
static_cast<float>(countingDetails.
m_nTotal) : 0.f)
521 <<
"%|, correctId " << ((countingDetails.
m_nTotal - countingDetails.
m_nMatch0 > 0) ? 100.f *
static_cast<float>(countingDetails.
m_correctId) /
static_cast<float>(countingDetails.
m_nTotal - countingDetails.
m_nMatch0) : 0.f)
536 std::ofstream mapFile, eventFile;
540 std::cout << std::endl <<
"EVENT INFO " << std::endl;
541 mapFile << std::endl <<
"EVENT INFO " << std::endl;
546 for (
const InteractionTargetResultMap::value_type &interactionMapEntry : interactionTargetResultMap)
551 unsigned int nCorrectEvents(0);
553 for (
const TargetResult &targetResult : targetResultList)
555 if (targetResult.m_isCorrect)
560 eventFile <<
"Correct event: fileId: " << targetResult.m_fileIdentifier <<
", eventNumber: " << targetResult.m_eventNumber <<
", interactionType " <<
ToString(interactionType) << std::endl;
565 for (
const PrimaryResultMap::value_type &primaryMapEntry : primaryResultMap)
594 std::cout <<
ToString(interactionType) << std::endl <<
"-nEvents " << targetResultList.size() <<
", nCorrect " << nCorrectEvents
595 <<
", fCorrect " << 100.f *
static_cast<float>(nCorrectEvents) /
static_cast<float>(targetResultList.size()) <<
"%" << std::endl;
599 mapFile <<
ToString(interactionType) << std::endl <<
"-nEvents " << targetResultList.size() <<
", nCorrect " << nCorrectEvents
600 <<
", fCorrect " << 100.f *
static_cast<float>(nCorrectEvents) /
static_cast<float>(targetResultList.size()) <<
"%" << std::endl;
617 targetHistogramCollection.
m_hVtxDeltaX =
new TH1F((histPrefix +
"VtxDeltaX").c_str(),
"", 40000, -2000., 2000.);
618 targetHistogramCollection.
m_hVtxDeltaX->GetXaxis()->SetRangeUser(-5., +5.);
619 targetHistogramCollection.
m_hVtxDeltaX->GetXaxis()->SetTitle(
"Vertex #DeltaX [cm]");
620 targetHistogramCollection.
m_hVtxDeltaX->GetYaxis()->SetTitle(
"Number of Events");
625 targetHistogramCollection.
m_hVtxDeltaY =
new TH1F((histPrefix +
"VtxDeltaY").c_str(),
"", 40000, -2000., 2000.);
626 targetHistogramCollection.
m_hVtxDeltaY->GetXaxis()->SetRangeUser(-5., +5.);
627 targetHistogramCollection.
m_hVtxDeltaY->GetXaxis()->SetTitle(
"Vertex #DeltaY [cm]");
628 targetHistogramCollection.
m_hVtxDeltaY->GetYaxis()->SetTitle(
"Number of Events");
633 targetHistogramCollection.
m_hVtxDeltaZ =
new TH1F((histPrefix +
"VtxDeltaZ").c_str(),
"", 40000, -2000., 2000.);
634 targetHistogramCollection.
m_hVtxDeltaZ->GetXaxis()->SetRangeUser(-5., +5.);
635 targetHistogramCollection.
m_hVtxDeltaZ->GetXaxis()->SetTitle(
"Vertex #DeltaZ [cm]");
636 targetHistogramCollection.
m_hVtxDeltaZ->GetYaxis()->SetTitle(
"Number of Events");
641 targetHistogramCollection.
m_hVtxDeltaR =
new TH1F((histPrefix +
"VtxDeltaR").c_str(),
"", 40000, -100., 1900.);
642 targetHistogramCollection.
m_hVtxDeltaR->GetXaxis()->SetRangeUser(0., +5.);
643 targetHistogramCollection.
m_hVtxDeltaR->GetXaxis()->SetTitle(
"Vertex #DeltaR [cm]");
644 targetHistogramCollection.
m_hVtxDeltaR->GetYaxis()->SetTitle(
"Number of Events");
657 const int nHitBins(35);
const int nHitBinEdges(nHitBins + 1);
658 float hitsBinning[nHitBinEdges];
659 for (
int n = 0; n < nHitBinEdges; ++n) hitsBinning[n] = std::pow(10., 1 +
static_cast<float>(n + 2) / 10.);
663 primaryHistogramCollection.
m_hHitsAll =
new TH1F((histPrefix +
"HitsAll").c_str(),
"", nHitBins, hitsBinning);
664 primaryHistogramCollection.
m_hHitsAll->GetXaxis()->SetRangeUser(1., +6000);
665 primaryHistogramCollection.
m_hHitsAll->GetXaxis()->SetTitle(
"Number of Hits");
666 primaryHistogramCollection.
m_hHitsAll->GetYaxis()->SetTitle(
"Number of Events");
671 primaryHistogramCollection.
m_hHitsEfficiency =
new TH1F((histPrefix +
"HitsEfficiency").c_str(),
"", nHitBins, hitsBinning);
672 primaryHistogramCollection.
m_hHitsEfficiency->GetXaxis()->SetRangeUser(1., +6000);
673 primaryHistogramCollection.
m_hHitsEfficiency->GetXaxis()->SetTitle(
"Number of Hits");
674 primaryHistogramCollection.
m_hHitsEfficiency->GetYaxis()->SetRangeUser(0., +1.01);
675 primaryHistogramCollection.
m_hHitsEfficiency->GetYaxis()->SetTitle(
"Reconstruction Efficiency");
678 const int nMomentumBins(26);
const int nMomentumBinEdges(nMomentumBins + 1);
679 float momentumBinning[nMomentumBinEdges] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.4, 1.6, 2.0, 2.4, 2.8, 3.4, 4., 5., 10., 15., 20., 30., 40., 50.};
683 primaryHistogramCollection.
m_hMomentumAll =
new TH1F((histPrefix +
"MomentumAll").c_str(),
"", nMomentumBins, momentumBinning);
684 primaryHistogramCollection.
m_hMomentumAll->GetXaxis()->SetRangeUser(0., +5.5);
685 primaryHistogramCollection.
m_hMomentumAll->GetXaxis()->SetTitle(
"True Momentum [GeV]");
686 primaryHistogramCollection.
m_hMomentumAll->GetYaxis()->SetTitle(
"Number of Events");
691 primaryHistogramCollection.
m_hMomentumEfficiency =
new TH1F((histPrefix +
"MomentumEfficiency").c_str(),
"", nMomentumBins, momentumBinning);
700 primaryHistogramCollection.
m_hCompleteness =
new TH1F((histPrefix +
"Completeness").c_str(),
"", 51, -0.01, 1.01);
701 primaryHistogramCollection.
m_hCompleteness->GetXaxis()->SetTitle(
"Completeness");
702 primaryHistogramCollection.
m_hCompleteness->GetYaxis()->SetRangeUser(0., +1.01);
703 primaryHistogramCollection.
m_hCompleteness->GetYaxis()->SetTitle(
"Fraction of Events");
706 if (!primaryHistogramCollection.
m_hPurity)
708 primaryHistogramCollection.
m_hPurity =
new TH1F((histPrefix +
"Purity").c_str(),
"", 51, -0.01, 1.01);
709 primaryHistogramCollection.
m_hPurity->GetXaxis()->SetTitle(
"Purity");
710 primaryHistogramCollection.
m_hPurity->GetYaxis()->SetRangeUser(0., +1.01);
711 primaryHistogramCollection.
m_hPurity->GetYaxis()->SetTitle(
"Fraction of Events");
731 for (InteractionPrimaryHistogramMap::const_iterator iter = interactionPrimaryHistogramMap.begin(), iterEnd = interactionPrimaryHistogramMap.end(); iter != iterEnd; ++iter)
736 for (PrimaryHistogramMap::const_iterator hIter = primaryHistogramMap.begin(), hIterEnd = primaryHistogramMap.end(); hIter != hIterEnd; ++hIter)
741 for (
int n = -1; n <= primaryHistogramCollection.
m_hHitsEfficiency->GetXaxis()->GetNbins(); ++n)
743 const float found = primaryHistogramCollection.
m_hHitsEfficiency->GetBinContent(n + 1);
744 const float all = primaryHistogramCollection.
m_hHitsAll->GetBinContent(n + 1);
745 const float efficiency = (all > 0.f) ? found / all : 0.f;
746 const float error = (all > found) ? std::sqrt(efficiency * (1. - efficiency) / all) : 0.f;
754 const float all = primaryHistogramCollection.
m_hMomentumAll->GetBinContent(n + 1);
755 const float efficiency = (all > 0.f) ? found / all : 0.f;
756 const float error = (all > found) ? std::sqrt(efficiency * (1. - efficiency) / all) : 0.f;
762 primaryHistogramCollection.
m_hPurity->Scale(1. /
static_cast<double>(primaryHistogramCollection.
m_hPurity->GetEntries()));
771 switch (expectedPrimary)
773 case MUON :
return "MUON";
775 case PROTON1 :
return "PROTON1";
776 case PROTON2 :
return "PROTON2";
777 case PROTON3 :
return "PROTON3";
778 case PROTON4 :
return "PROTON4";
779 case PROTON5 :
return "PROTON5";
780 case PIPLUS :
return "PIPLUS";
781 case PIMINUS :
return "PIMINUS";
782 case NEUTRON :
return "NEUTRON";
783 case PHOTON1 :
return "PHOTON1";
784 case PHOTON2 :
return "PHOTON2";
786 default:
return "UNKNOWN";
794 switch (interactionType)
802 case CCQEL_E:
return "CCQEL_E";
808 case NCQEL_P:
return "NCQEL_P";
837 case CCRES_E:
return "CCRES_E";
861 case NCRES_P:
return "NCRES_P";
914 case NCDIS_P:
return "NCDIS_P";
943 case CCCOH:
return "CCCOH";
944 case NCCOH:
return "NCCOH";
961 default:
return "UNKNOWN";
bool PassSBNDFiducialCut(const SimpleMCTarget &simpleMCTarget)
Whether a simple mc event passes sbnd fiducial cut, applied to target vertices.
bool IsGoodParticleIdMatch(const SimpleMCPrimary &simpleMCPrimary, const int bestMatchPfoPdgCode)
Whether a provided mc primary and best matched pfo are deemed to have a good particle id match.
void DisplayInteractionCountingMap(const InteractionCountingMap &interactionCountingMap, const Parameters ¶meters)
Print details to screen for a provided interaction type to counting map.
void CountPfoMatches(const SimpleMCEvent &simpleMCEvent, const Parameters ¶meters, InteractionCountingMap &interactionCountingMap, InteractionTargetResultMap &interactionTargetResultMap)
CountPfoMatches Relies on fact that primary list is sorted by number of true good hits.
int ReadNextEvent(TChain *const pTChain, const int iEntry, SimpleMCEvent &simpleMCEvent, const Parameters ¶meters)
Read the next event from the chain.
ExpectedPrimary GetExpectedPrimary(const SimpleMCPrimary &simpleMCPrimary, const SimpleMCPrimaryList &simpleMCPrimaryList)
Work out which of the primary particles (expected for a given interaction types) corresponds to the p...
bool PassUbooneFiducialCut(const SimpleMCTarget &simpleMCTarget)
Whether a simple mc event passes uboone fiducial cut, applied to target vertices.
bool PassFiducialCut(const SimpleMCTarget &simpleMCTarget, const Parameters ¶meters)
Whether a simple mc event passes the relevant fiducial cut, applied to target vertices.
void DisplaySimpleMCEventMatches(const SimpleMCEvent &simpleMCEvent, const Parameters ¶meters)
Print matching details to screen for a simple mc event.
void FillTargetHistogramCollection(const std::string &histPrefix, const TargetResult &targetResult, TargetHistogramCollection &targetHistogramCollection)
Fill histograms in the provided target histogram collection, using information in the provided target...
void ProcessHistogramCollections(const InteractionPrimaryHistogramMap &interactionPrimaryHistogramMap)
Process histograms stored in the provided map e.g. calculating final efficiencies,...
void FillPrimaryHistogramCollection(const std::string &histPrefix, const Parameters ¶meters, const PrimaryResult &primaryResult, PrimaryHistogramCollection &primaryHistogramCollection)
Fill histograms in the provided histogram collection, using information in the provided primary resul...
void AnalyseInteractionTargetResultMap(const InteractionTargetResultMap &interactionTargetResultMap, const Parameters ¶meters)
Opportunity to fill histograms, perform post-processing of information collected in main loop over nt...
void Validation(const std::string &inputFiles, const Parameters ¶meters)
Validation - Main entry point for analysis.
std::string ToString(const ExpectedPrimary expectedPrimary)
Get a string representation of an interaction type.
std::map< InteractionType, PrimaryHistogramMap > InteractionPrimaryHistogramMap
std::map< ExpectedPrimary, CountingDetails > CountingMap
std::map< ExpectedPrimary, PrimaryHistogramCollection > PrimaryHistogramMap
std::map< InteractionType, CountingMap > InteractionCountingMap
std::map< ExpectedPrimary, PrimaryResult > PrimaryResultMap
std::map< InteractionType, TargetResultList > InteractionTargetResultMap
std::vector< TargetResult > TargetResultList
std::vector< float > FloatVector
InteractionType
InteractionType enum.
@ CCRES_E_P_P_P_P_P_PIZERO
@ BEAM_PARTICLE_KAON_MINUS
@ CCRES_E_P_P_P_P_P_PIPLUS
@ CCRES_MU_P_P_P_P_PIPLUS
@ CCDIS_MU_P_P_P_P_PIPLUS
@ BEAM_PARTICLE_KAON_PLUS
@ CCRES_MU_P_P_P_P_P_PHOTON
@ CCDIS_MU_P_P_P_P_P_PIPLUS
@ CCRES_MU_P_P_P_P_P_PIPLUS
@ CCDIS_MU_P_P_P_P_PHOTON
@ CCRES_E_P_P_P_P_P_PHOTON
@ CCRES_MU_P_P_P_P_P_PIZERO
@ CCRES_MU_P_P_P_P_PHOTON
@ NCDIS_P_P_P_P_P_PIMINUS
@ CCDIS_MU_P_P_P_P_P_PIZERO
@ CCRES_MU_P_P_P_P_PIZERO
@ CCDIS_MU_P_P_P_P_PIZERO
@ NCRES_P_P_P_P_P_PIMINUS
@ CCDIS_MU_P_P_P_P_P_PHOTON
std::vector< int > IntVector
ExpectedPrimary
ExpectedPrimary enum.
std::vector< SimpleMCPrimary > SimpleMCPrimaryList
std::map< InteractionType, TargetHistogramCollection > InteractionTargetHistogramMap
unsigned int m_nTotal
The total number of occurences.
unsigned int m_correctId
The number of times the mc primary particle id was correct.
unsigned int m_nMatch0
The number of times the mc primary has 0 pfo matches.
unsigned int m_nMatch2
The number of times the mc primary has 2 pfo matches.
unsigned int m_nMatch3Plus
The number of times the mc primary has 3 or more pfo matches.
unsigned int m_nMatch1
The number of times the mc primary has 1 pfo matches.
bool m_applyUbooneFiducialCut
Whether to apply uboone fiducial volume cut to true neutrino vertex position.
bool m_triggeredBeamOnly
Whether to only consider triggered beam particles.
std::string m_eventFileName
File name to which to write list of correct events.
std::string m_histPrefix
Histogram name prefix.
float m_vertexXCorrection
The vertex x correction, added to reported mc neutrino endpoint x value, in cm.
std::string m_mapFileName
File name to which to write output ascii tables, etc.
int m_nEventsToProcess
The number of events to process.
bool m_displayMatchedEvents
Whether to display matching results for individual events.
bool m_testBeamMode
Whether running in test beam mode.
int m_skipEvents
The number of events to skip.
bool m_histogramOutput
Whether to produce output histograms.
bool m_correctTrackShowerId
Whether to demand that pfos are correctly flagged as tracks or showers.
bool m_applySBNDFiducialCut
Whether to apply sbnd fiducial volume cut to true neutrino vertex position.
PrimaryHistogramCollection class.
TH1F * m_hMomentumAll
The number of primaries vs momentum histogram.
TH1F * m_hHitsEfficiency
The primary efficiency vs number of hits histogram.
TH1F * m_hPurity
The primary (best match) purity histogram.
TH1F * m_hMomentumEfficiency
The primary efficiency vs momentum histogram.
TH1F * m_hHitsAll
The number of primaries vs number of hits histogram.
TH1F * m_hCompleteness
The primary (best match) completeness histogram.
bool m_isCorrectParticleId
Whether the best matched pfo has the correct particle id.
unsigned int m_nPfoMatches
The total number of pfo matches for a given primary.
unsigned int m_nMCHitsTotal
The number of hits in the mc primary.
unsigned int m_nBestMatchRecoHitsTotal
The number of hits in the best matched pfo.
float m_bestMatchPurity
The purity of the best matched pfo.
float m_trueMomentum
The true momentum of the mc primary.
float m_bestMatchCompleteness
The completeness of the best matched pfo.
unsigned int m_nBestMatchSharedHitsTotal
The number of hits shared by the mc primary and the best matched pfo.
int m_eventNumber
The event number.
int m_fileIdentifier
The file identifier.
SimpleMCTargetList m_mcTargetList
The list of mc targets.
int m_nMCTargets
The number of mc targets.
SimpleThreeVector m_momentum
The momentum.
int m_nMCHitsV
The number of v mc hits.
int m_pdgCode
The pdg code.
int m_bestMatchPfoIsTestBeam
Whether best match pfo is reconstructed as a test beam particle.
int m_bestMatchPfoNHitsW
The best match pfo number of w pfo hits.
int m_primaryId
The identifier.
int m_bestMatchPfoNSharedHitsW
The best match pfo number of w matched hits.
int m_bestMatchPfoIsRecoNu
Whether best match pfo is reconstructed as part of a neutrino hierarchy.
int m_bestMatchPfoRecoNuId
The identifier of the associated reco neutrino (if part of a neutrino hierarchy)
int m_nPrimaryMatchedPfos
The number of matched pfos.
int m_bestMatchPfoNSharedHitsU
The best match pfo number of u matched hits.
int m_bestMatchPfoNHitsTotal
The best match pfo total number of pfo hits.
int m_nMCHitsTotal
The total number of mc hits.
SimpleThreeVector m_vertex
The vertex.
int m_bestMatchPfoNHitsU
The best match pfo number of u pfo hits.
int m_nMCHitsU
The number of u mc hits.
int m_bestMatchPfoNHitsV
The best match pfo number of v pfo hits.
SimpleThreeVector m_endpoint
The endpoint.
int m_bestMatchPfoId
The best match pfo identifier.
int m_bestMatchPfoNSharedHitsV
The best match pfo number of v matched hits.
int m_nPrimaryMatchedCRPfos
The number of matched cr pfos.
float m_energy
The energy.
int m_nMCHitsW
The number of w mc hits.
int m_bestMatchPfoPdgCode
The best match pfo pdg code.
int m_nPrimaryMatchedNuPfos
The number of matched nu pfos.
int m_bestMatchPfoNSharedHitsTotal
The best match pfo total number of matched hits.
int m_isCosmicRay
Whether the target is a cosmic ray.
int m_isFakeCR
Whether the target was reconstructed as a fake cosmic ray.
int m_isNeutrino
Whether the target is a neutrino.
int m_nTargetNuMatches
The number of neutrino pfo matches to the target.
int m_isLost
Whether the target was lost (not reconstructed)
int m_nTargetGoodNuMatches
The number of good neutrino pfo matches to the target (all from same parent neutrino)
int m_isFakeNu
Whether the target was reconstructed as a fake neutrino.
SimpleMCPrimaryList m_mcPrimaryList
The list of mc primaries.
int m_nTargetNuLosses
The number of neutrino primaries with no matches.
int m_isCorrectNu
Whether the target was correctly reconstructed as a neutrino.
int m_interactionType
The target interaction type.
int m_nTargetPrimaries
The number of target mc primaries.
int m_isCorrectTB
Whether the target was correctly reconstructed as a beam particle.
int m_nTargetMatches
The number of pfo matches to the target.
int m_nTargetCRMatches
The number of cosmic ray pfo matches to the target.
SimpleThreeVector m_targetVertex
The target vertex position.
int m_nTargetNuSplits
The number of split neutrino pfo matches to the target (from different parent neutrinos)
int m_isSplitNu
Whether the target was reconstructed as a split neutrino.
int m_isCorrectCR
Whether the target was correctly reconstructed as a cosmic ray.
SimpleThreeVector m_recoVertex
The reco vertex position, if available.
int m_mcNuanceCode
The target nuance code.
int m_isBeamParticle
Whether the target is a beam particle.
int m_isSplitCR
Whether the target was reconstructed as a split cosmic ray.
TargetHistogramCollection class.
TH1F * m_hVtxDeltaY
The vtx delta y histogram.
TH1F * m_hVtxDeltaX
The vtx delta x histogram.
TH1F * m_hVtxDeltaR
The vtx delta r histogram.
TH1F * m_hVtxDeltaZ
The vtx delta z histogram.
bool m_hasRecoVertex
Whether a reco vertex is matched to the target.
bool m_isCorrect
Whether the target is reconstructed correctly.
PrimaryResultMap m_primaryResultMap
The primary result map.
int m_fileIdentifier
The file identifier.
int m_eventNumber
The event number.
SimpleThreeVector m_vertexOffset
The offset between the reco and true target vertices.