Pandora
Pandora source code navigator
Loading...
Searching...
No Matches
Validation.C
Go to the documentation of this file.
1
8#include "TChain.h"
9#include "TH1F.h"
10
11#include "Validation.h"
12
13#include <iomanip>
14#include <iostream>
15#include <fstream>
16#include <sstream>
17
18void Validation(const std::string &inputFiles, const Parameters &parameters)
19{
20 TChain *pTChain = new TChain("Validation", "pTChain");
21 pTChain->Add(inputFiles.c_str());
22
23 InteractionCountingMap interactionCountingMap;
24 InteractionTargetResultMap interactionTargetResultMap;
25
26 int nEvents(0), nProcessedEvents(0);
27 const int nChainEntries(pTChain->GetEntries());
28
29 for (int iEntry = 0; iEntry < nChainEntries; )
30 {
31 SimpleMCEvent simpleMCEvent;
32 iEntry += ReadNextEvent(pTChain, iEntry, simpleMCEvent, parameters);
33
34 if (nEvents++ < parameters.m_skipEvents)
35 continue;
36
37 if (nEvents % 50 == 0)
38 std::cout << "nEvents " << nEvents << "\r" << std::flush;
39
40 if (nProcessedEvents++ >= parameters.m_nEventsToProcess)
41 break;
42
43 if (parameters.m_displayMatchedEvents)
44 DisplaySimpleMCEventMatches(simpleMCEvent, parameters);
45
46 CountPfoMatches(simpleMCEvent, parameters, interactionCountingMap, interactionTargetResultMap);
47 }
48
49 DisplayInteractionCountingMap(interactionCountingMap, parameters);
50 AnalyseInteractionTargetResultMap(interactionTargetResultMap, parameters);
51}
52
53//------------------------------------------------------------------------------------------------------------------------------------------
54
55int ReadNextEvent(TChain *const pTChain, const int iEntry, SimpleMCEvent &simpleMCEvent, const Parameters &parameters)
56{
57 int thisEventNumber(0), iTarget(0);
58 const int nChainEntries(pTChain->GetEntries());
59
60 pTChain->SetBranchAddress("eventNumber", &thisEventNumber);
61 pTChain->SetBranchAddress("fileIdentifier", &simpleMCEvent.m_fileIdentifier);
62 pTChain->GetEntry(iEntry);
63 simpleMCEvent.m_eventNumber = thisEventNumber;
64
65 while (iEntry + iTarget < nChainEntries)
66 {
67 SimpleMCTarget simpleMCTarget;
68
69 pTChain->SetBranchAddress("interactionType", &simpleMCTarget.m_interactionType);
70 pTChain->SetBranchAddress("mcNuanceCode", &simpleMCTarget.m_mcNuanceCode);
71 pTChain->SetBranchAddress("isCosmicRay", &simpleMCTarget.m_isCosmicRay);
72 pTChain->SetBranchAddress("targetVertexX", &simpleMCTarget.m_targetVertex.m_x);
73 pTChain->SetBranchAddress("targetVertexY", &simpleMCTarget.m_targetVertex.m_y);
74 pTChain->SetBranchAddress("targetVertexZ", &simpleMCTarget.m_targetVertex.m_z);
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);
83 pTChain->SetBranchAddress("nTargetCRMatches", &simpleMCTarget.m_nTargetCRMatches);
84 pTChain->SetBranchAddress("nTargetPrimaries", &simpleMCTarget.m_nTargetPrimaries);
85
86 if (parameters.m_testBeamMode)
87 {
88 pTChain->SetBranchAddress("isBeamParticle", &simpleMCTarget.m_isBeamParticle);
89 pTChain->SetBranchAddress("isCorrectTB", &simpleMCTarget.m_isCorrectTB);
90 }
91 else
92 {
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);
97 pTChain->SetBranchAddress("nTargetNuMatches", &simpleMCTarget.m_nTargetNuMatches);
98 }
99
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);
107
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);
136
137 if (parameters.m_testBeamMode)
138 {
139 pTChain->SetBranchAddress("bestMatchPfoIsTB", &pBestMatchPfoIsTestBeam);
140 }
141 else
142 {
143 pTChain->SetBranchAddress("nPrimaryMatchedNuPfos", &pNPrimaryMatchedNuPfos);
144 pTChain->SetBranchAddress("bestMatchPfoIsRecoNu", &pBestMatchPfoIsRecoNu);
145 pTChain->SetBranchAddress("bestMatchPfoRecoNuId", &pBestMatchPfoRecoNuId);
146 pTChain->SetBranchAddress("nTargetGoodNuMatches", &simpleMCTarget.m_nTargetGoodNuMatches);
147 pTChain->SetBranchAddress("nTargetNuSplits", &simpleMCTarget.m_nTargetNuSplits);
148 pTChain->SetBranchAddress("nTargetNuLosses", &simpleMCTarget.m_nTargetNuLosses);
149 }
150
151 pTChain->GetEntry(iEntry + iTarget++);
152
153 if (simpleMCEvent.m_eventNumber != thisEventNumber)
154 break;
155
156 for (int iPrimary = 0; iPrimary < simpleMCTarget.m_nTargetPrimaries; ++iPrimary)
157 {
158 SimpleMCPrimary simpleMCPrimary;
159 simpleMCPrimary.m_primaryId = pMCPrimaryId->at(iPrimary);
160 simpleMCPrimary.m_pdgCode = pMCPrimaryPdg->at(iPrimary);
161 simpleMCPrimary.m_energy = pMCPrimaryE->at(iPrimary);
162 simpleMCPrimary.m_momentum.m_x = pMCPrimaryPX->at(iPrimary);
163 simpleMCPrimary.m_momentum.m_y = pMCPrimaryPY->at(iPrimary);
164 simpleMCPrimary.m_momentum.m_z = pMCPrimaryPZ->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);
171 simpleMCPrimary.m_nMCHitsTotal = pNMCHitsTotal->at(iPrimary);
172 simpleMCPrimary.m_nMCHitsU = pNMCHitsU->at(iPrimary);
173 simpleMCPrimary.m_nMCHitsV = pNMCHitsV->at(iPrimary);
174 simpleMCPrimary.m_nMCHitsW = pNMCHitsW->at(iPrimary);
175 simpleMCPrimary.m_nPrimaryMatchedPfos = pNPrimaryMatchedPfos->at(iPrimary);
176 simpleMCPrimary.m_nPrimaryMatchedCRPfos = pNPrimaryMatchedCRPfos->at(iPrimary);
177 simpleMCPrimary.m_bestMatchPfoId = pBestMatchPfoId->at(iPrimary);
178 simpleMCPrimary.m_bestMatchPfoPdgCode = pBestMatchPfoPdg->at(iPrimary);
179 simpleMCPrimary.m_bestMatchPfoNHitsTotal = pBestMatchPfoNHitsTotal->at(iPrimary);
180 simpleMCPrimary.m_bestMatchPfoNHitsU = pBestMatchPfoNHitsU->at(iPrimary);
181 simpleMCPrimary.m_bestMatchPfoNHitsV = pBestMatchPfoNHitsV->at(iPrimary);
182 simpleMCPrimary.m_bestMatchPfoNHitsW = pBestMatchPfoNHitsW->at(iPrimary);
183 simpleMCPrimary.m_bestMatchPfoNSharedHitsTotal = pBestMatchPfoNSharedHitsTotal->at(iPrimary);
184 simpleMCPrimary.m_bestMatchPfoNSharedHitsU = pBestMatchPfoNSharedHitsU->at(iPrimary);
185 simpleMCPrimary.m_bestMatchPfoNSharedHitsV = pBestMatchPfoNSharedHitsV->at(iPrimary);
186 simpleMCPrimary.m_bestMatchPfoNSharedHitsW = pBestMatchPfoNSharedHitsW->at(iPrimary);
187
188 if (parameters.m_testBeamMode)
189 {
190 simpleMCPrimary.m_bestMatchPfoIsTestBeam = pBestMatchPfoIsTestBeam->at(iPrimary);
191 }
192 else
193 {
194 simpleMCPrimary.m_nPrimaryMatchedNuPfos = pNPrimaryMatchedNuPfos->at(iPrimary);
195 simpleMCPrimary.m_bestMatchPfoIsRecoNu = pBestMatchPfoIsRecoNu->at(iPrimary);
196 simpleMCPrimary.m_bestMatchPfoRecoNuId = pBestMatchPfoRecoNuId->at(iPrimary);
197 }
198
199 simpleMCTarget.m_mcPrimaryList.push_back(simpleMCPrimary);
200 }
201
202 simpleMCEvent.m_mcTargetList.push_back(simpleMCTarget);
203 simpleMCEvent.m_nMCTargets = simpleMCEvent.m_mcTargetList.size();
204 }
205
206 pTChain->ResetBranchAddresses();
207 return simpleMCEvent.m_nMCTargets;
208}
209
210//------------------------------------------------------------------------------------------------------------------------------------------
211
212void DisplaySimpleMCEventMatches(const SimpleMCEvent &simpleMCEvent, const Parameters &parameters)
213{
214 std::cout << "---INTERPRETED-MATCHING-OUTPUT------------------------------------------------------------------" << std::endl;
215 std::cout << "File " << simpleMCEvent.m_fileIdentifier << ", event " << simpleMCEvent.m_eventNumber << std::endl;
216
217 int nCorrectNu(0), nTotalNu(0), nCorrectTB(0), nTotalTB(0), nCorrectCR(0), nTotalCR(0), nFakeNu(0), nFakeCR(0), nSplitNu(0), nSplitCR(0), nLost(0);
218
219 for (const SimpleMCTarget &simpleMCTarget : simpleMCEvent.m_mcTargetList)
220 {
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;
225
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;
240
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;
252
253 for (const SimpleMCPrimary &simpleMCPrimary : simpleMCTarget.m_mcPrimaryList)
254 {
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;
269
270 if (0 == simpleMCPrimary.m_nPrimaryMatchedPfos)
271 {
272 std::cout << "-No matched Pfo" << std::endl;
273 continue;
274 }
275
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;
291 }
292 }
293
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;
307}
308
309//------------------------------------------------------------------------------------------------------------------------------------------
310
311void CountPfoMatches(const SimpleMCEvent &simpleMCEvent, const Parameters &parameters, InteractionCountingMap &interactionCountingMap,
312 InteractionTargetResultMap &interactionTargetResultMap)
313{
314 for (const SimpleMCTarget &simpleMCTarget : simpleMCEvent.m_mcTargetList)
315 {
316 if ((!PassFiducialCut(simpleMCTarget, parameters) && simpleMCTarget.m_isNeutrino) ||
317 (parameters.m_triggeredBeamOnly && simpleMCTarget.m_isBeamParticle && simpleMCTarget.m_mcNuanceCode != 2001))
318 continue;
319
320 TargetResult targetResult;
321 targetResult.m_fileIdentifier = simpleMCEvent.m_fileIdentifier;
322 targetResult.m_eventNumber = simpleMCEvent.m_eventNumber;
323 targetResult.m_isCorrect = (simpleMCTarget.m_isNeutrino && simpleMCTarget.m_isCorrectNu) ||
324 (simpleMCTarget.m_isBeamParticle && simpleMCTarget.m_isCorrectTB) ||
325 (simpleMCTarget.m_isCosmicRay && simpleMCTarget.m_isCorrectCR);
326
327 if (simpleMCTarget.m_nTargetMatches > 0)
328 {
329 targetResult.m_hasRecoVertex = true;
330 targetResult.m_vertexOffset = simpleMCTarget.m_recoVertex - simpleMCTarget.m_targetVertex;
331 targetResult.m_vertexOffset.m_x = targetResult.m_vertexOffset.m_x - parameters.m_vertexXCorrection;
332 }
333
334 const InteractionType interactionType(static_cast<InteractionType>(simpleMCTarget.m_interactionType));
335
336 for (const SimpleMCPrimary &simpleMCPrimary : simpleMCTarget.m_mcPrimaryList)
337 {
338 const ExpectedPrimary expectedPrimary(GetExpectedPrimary(simpleMCPrimary, simpleMCTarget.m_mcPrimaryList));
339
340 PrimaryResult &primaryResult = targetResult.m_primaryResultMap[expectedPrimary];
341 CountingDetails &countingDetails = interactionCountingMap[interactionType][expectedPrimary];
342 ++countingDetails.m_nTotal;
343
344 // ATTN Fail cosmic ray matches to neutrinos (or beam particles) and vice versa
345 bool incorrectMatchToCR(parameters.m_testBeamMode ? (simpleMCTarget.m_isCosmicRay == simpleMCPrimary.m_bestMatchPfoIsTestBeam) : (simpleMCTarget.m_isCosmicRay == simpleMCPrimary.m_bestMatchPfoIsRecoNu));
346
347 if ((simpleMCPrimary.m_bestMatchPfoId >= 0) && incorrectMatchToCR)
348 {
349 ++countingDetails.m_nMatch0;
350 continue;
351 }
352
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;
356 else ++countingDetails.m_nMatch3Plus;
357
358 primaryResult.m_nPfoMatches = simpleMCPrimary.m_nPrimaryMatchedPfos;
359 primaryResult.m_nMCHitsTotal = simpleMCPrimary.m_nMCHitsTotal;
360 primaryResult.m_nBestMatchSharedHitsTotal = simpleMCPrimary.m_bestMatchPfoNSharedHitsTotal;
361 primaryResult.m_nBestMatchRecoHitsTotal = simpleMCPrimary.m_bestMatchPfoNHitsTotal;
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;
364 primaryResult.m_isCorrectParticleId = IsGoodParticleIdMatch(simpleMCPrimary, simpleMCPrimary.m_bestMatchPfoPdgCode);
365
366 if ((simpleMCPrimary.m_nPrimaryMatchedPfos > 0) && primaryResult.m_isCorrectParticleId)
367 ++countingDetails.m_correctId;
368
369 if (parameters.m_correctTrackShowerId && !primaryResult.m_isCorrectParticleId)
370 targetResult.m_isCorrect = false;
371
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));
373 primaryResult.m_trueMomentum = pTot;
374 }
375
376 interactionTargetResultMap[interactionType].push_back(targetResult);
377 }
378}
379
380//------------------------------------------------------------------------------------------------------------------------------------------
381
382bool PassFiducialCut(const SimpleMCTarget &simpleMCTarget, const Parameters &parameters)
383{
384 if (parameters.m_applyUbooneFiducialCut && parameters.m_applySBNDFiducialCut)
385 throw std::invalid_argument("Parameters has fiducial cuts for uBooNE and SBND");
386
387 if (parameters.m_applyUbooneFiducialCut)
388 return PassUbooneFiducialCut(simpleMCTarget);
389
390 if (parameters.m_applySBNDFiducialCut)
391 return PassSBNDFiducialCut(simpleMCTarget);
392
393 return true;
394}
395
396//------------------------------------------------------------------------------------------------------------------------------------------
397
398bool PassUbooneFiducialCut(const SimpleMCTarget &simpleMCTarget)
399{
400 const float eVx(256.35), eVy(233.), eVz(1036.8);
401 const float xBorder(10.), yBorder(20.), zBorder(10.);
402
403 if ((simpleMCTarget.m_targetVertex.m_x < (eVx - xBorder)) && (simpleMCTarget.m_targetVertex.m_x > xBorder) &&
404 (simpleMCTarget.m_targetVertex.m_y < (eVy / 2. - yBorder)) && (simpleMCTarget.m_targetVertex.m_y > (-eVy / 2. + yBorder)) &&
405 (simpleMCTarget.m_targetVertex.m_z < (eVz - zBorder)) && (simpleMCTarget.m_targetVertex.m_z > zBorder) )
406 {
407 return true;
408 }
409
410 return false;
411}
412
413//------------------------------------------------------------------------------------------------------------------------------------------
414
415bool PassSBNDFiducialCut(const SimpleMCTarget &simpleMCTarget)
416{
417 const float eVx(400.f), eVy(400.f), eVz(500.f);
418 const float xBorder(10.f), yBorder(20.f), zBorder(10.f);
419
420 // ATTN origin definition is different in SBND to uBooNE. Both x & y are centered in the middle of the face
421 if ((simpleMCTarget.m_targetVertex.m_x < (eVx / 2. - xBorder)) && (simpleMCTarget.m_targetVertex.m_x > (-eVx / 2. + xBorder)) &&
422 (simpleMCTarget.m_targetVertex.m_y < (eVy / 2. - yBorder)) && (simpleMCTarget.m_targetVertex.m_y > (-eVy / 2. + yBorder)) &&
423 (simpleMCTarget.m_targetVertex.m_z < (eVz - zBorder)) && (simpleMCTarget.m_targetVertex.m_z > zBorder) )
424 {
425 return true;
426 }
427
428 return false;
429}
430
431//------------------------------------------------------------------------------------------------------------------------------------------
432
433ExpectedPrimary GetExpectedPrimary(const SimpleMCPrimary &simpleMCPrimary, const SimpleMCPrimaryList &simpleMCPrimaryList)
434{
435 // ATTN: Relies on fact that primary list is sorted by number of good true hits
436 unsigned int nMuons(0), nElectrons(0), nProtons(0), nPiPlus(0), nPiMinus(0), nNeutrons(0), nPhotons(0);
437
438 for (const SimpleMCPrimary &simpleMCPrimaryInList : simpleMCPrimaryList)
439 {
440 if (&simpleMCPrimary == &simpleMCPrimaryInList)
441 {
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;
453 }
454
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;
462 }
463
464 return OTHER_PRIMARY;
465}
466
467//------------------------------------------------------------------------------------------------------------------------------------------
468
469bool IsGoodParticleIdMatch(const SimpleMCPrimary &simpleMCPrimary, const int bestMatchPfoPdgCode)
470{
471 const unsigned int absMCPdgCode(std::fabs(simpleMCPrimary.m_pdgCode));
472
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))) )
475 {
476 return false;
477 }
478
479 return true;
480}
481
482//------------------------------------------------------------------------------------------------------------------------------------------
483
484void DisplayInteractionCountingMap(const InteractionCountingMap &interactionCountingMap, const Parameters &parameters)
485{
486 std::cout << std::fixed;
487 std::cout << std::setprecision(1);
488
489 std::ofstream mapFile;
490 if (!parameters.m_mapFileName.empty()) mapFile.open(parameters.m_mapFileName, ios::app);
491
492 for (const InteractionCountingMap::value_type &interactionTypeMapEntry : interactionCountingMap)
493 {
494 const InteractionType interactionType(interactionTypeMapEntry.first);
495 const CountingMap &countingMap(interactionTypeMapEntry.second);
496 std::cout << std::endl << ToString(interactionType) << std::endl;
497
498 if (!parameters.m_mapFileName.empty())
499 mapFile << std::endl << ToString(interactionType) << std::endl;
500
501 for (const CountingMap::value_type &countingMapEntry : countingMap)
502 {
503 const ExpectedPrimary expectedPrimary(countingMapEntry.first);
504 const CountingDetails &countingDetails(countingMapEntry.second);
505
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)
512 << "%" << std::endl;
513
514 if (!parameters.m_mapFileName.empty())
515 {
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)
522 << "%" << std::endl;
523 }
524 }
525 }
526
527 if (!parameters.m_mapFileName.empty())
528 mapFile.close();
529}
530
531//------------------------------------------------------------------------------------------------------------------------------------------
532
533void AnalyseInteractionTargetResultMap(const InteractionTargetResultMap &interactionTargetResultMap, const Parameters &parameters)
534{
535 // Intended for filling histograms, post-processing of information collected in main loop over ntuple, etc.
536 std::ofstream mapFile, eventFile;
537 if (!parameters.m_mapFileName.empty()) mapFile.open(parameters.m_mapFileName, ios::app);
538 if (!parameters.m_eventFileName.empty()) eventFile.open(parameters.m_eventFileName, ios::app);
539
540 std::cout << std::endl << "EVENT INFO " << std::endl;
541 mapFile << std::endl << "EVENT INFO " << std::endl;
542
543 InteractionPrimaryHistogramMap interactionPrimaryHistogramMap;
544 InteractionTargetHistogramMap interactionTargetHistogramMap;
545
546 for (const InteractionTargetResultMap::value_type &interactionMapEntry : interactionTargetResultMap)
547 {
548 const InteractionType interactionType(interactionMapEntry.first);
549 const TargetResultList &targetResultList(interactionMapEntry.second);
550
551 unsigned int nCorrectEvents(0);
552
553 for (const TargetResult &targetResult : targetResultList)
554 {
555 if (targetResult.m_isCorrect)
556 {
557 ++nCorrectEvents;
558
559 if (!parameters.m_eventFileName.empty())
560 eventFile << "Correct event: fileId: " << targetResult.m_fileIdentifier << ", eventNumber: " << targetResult.m_eventNumber << ", interactionType " << ToString(interactionType) << std::endl;
561 }
562
563 const PrimaryResultMap &primaryResultMap(targetResult.m_primaryResultMap);
564
565 for (const PrimaryResultMap::value_type &primaryMapEntry : primaryResultMap)
566 {
567 const ExpectedPrimary expectedPrimary(primaryMapEntry.first);
568 const PrimaryResult &primaryResult(primaryMapEntry.second);
569
570 if (parameters.m_histogramOutput)
571 {
572 const std::string histPrefix(parameters.m_histPrefix + ToString(interactionType) + "_" + ToString(expectedPrimary) + "_");
573 PrimaryHistogramCollection &histogramCollection(interactionPrimaryHistogramMap[interactionType][expectedPrimary]);
574 FillPrimaryHistogramCollection(histPrefix, parameters, primaryResult, histogramCollection);
575
576 const std::string histPrefixAll(parameters.m_histPrefix + ToString(ALL_INTERACTIONS) + "_" + ToString(expectedPrimary) + "_");
577 PrimaryHistogramCollection &histogramCollectionAll(interactionPrimaryHistogramMap[ALL_INTERACTIONS][expectedPrimary]);
578 FillPrimaryHistogramCollection(histPrefixAll, parameters, primaryResult, histogramCollectionAll);
579 }
580 }
581
582 if (parameters.m_histogramOutput)
583 {
584 const std::string histPrefix(parameters.m_histPrefix + ToString(interactionType) + "_");
585 TargetHistogramCollection &histogramCollection(interactionTargetHistogramMap[interactionType]);
586 FillTargetHistogramCollection(histPrefix, targetResult, histogramCollection);
587
588 const std::string histPrefixAll(parameters.m_histPrefix + ToString(ALL_INTERACTIONS) + "_");
589 TargetHistogramCollection &histogramCollectionAll(interactionTargetHistogramMap[ALL_INTERACTIONS]);
590 FillTargetHistogramCollection(histPrefixAll, targetResult, histogramCollectionAll);
591 }
592 }
593
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;
596
597 if (!parameters.m_mapFileName.empty())
598 {
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;
601 }
602 }
603
604 if (parameters.m_histogramOutput)
605 ProcessHistogramCollections(interactionPrimaryHistogramMap);
606
607 if (!parameters.m_mapFileName.empty()) mapFile.close();
608 if (!parameters.m_eventFileName.empty()) eventFile.close();
609}
610
611//------------------------------------------------------------------------------------------------------------------------------------------
612
613void FillTargetHistogramCollection(const std::string &histPrefix, const TargetResult &targetResult, TargetHistogramCollection &targetHistogramCollection)
614{
615 if (!targetHistogramCollection.m_hVtxDeltaX)
616 {
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");
621 }
622
623 if (!targetHistogramCollection.m_hVtxDeltaY)
624 {
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");
629 }
630
631 if (!targetHistogramCollection.m_hVtxDeltaZ)
632 {
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");
637 }
638
639 if (!targetHistogramCollection.m_hVtxDeltaR)
640 {
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");
645 }
646
647 targetHistogramCollection.m_hVtxDeltaX->Fill(targetResult.m_vertexOffset.m_x);
648 targetHistogramCollection.m_hVtxDeltaY->Fill(targetResult.m_vertexOffset.m_y);
649 targetHistogramCollection.m_hVtxDeltaZ->Fill(targetResult.m_vertexOffset.m_z);
650 targetHistogramCollection.m_hVtxDeltaR->Fill(std::sqrt(targetResult.m_vertexOffset.m_x * targetResult.m_vertexOffset.m_x + targetResult.m_vertexOffset.m_y * targetResult.m_vertexOffset.m_y + targetResult.m_vertexOffset.m_z * targetResult.m_vertexOffset.m_z));
651}
652
653//------------------------------------------------------------------------------------------------------------------------------------------
654
655void FillPrimaryHistogramCollection(const std::string &histPrefix, const Parameters &parameters, const PrimaryResult &primaryResult, PrimaryHistogramCollection &primaryHistogramCollection)
656{
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.);
660
661 if (!primaryHistogramCollection.m_hHitsAll)
662 {
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");
667 }
668
669 if (!primaryHistogramCollection.m_hHitsEfficiency)
670 {
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");
676 }
677
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.};
680
681 if (!primaryHistogramCollection.m_hMomentumAll)
682 {
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");
687 }
688
689 if (!primaryHistogramCollection.m_hMomentumEfficiency)
690 {
691 primaryHistogramCollection.m_hMomentumEfficiency = new TH1F((histPrefix + "MomentumEfficiency").c_str(), "", nMomentumBins, momentumBinning);
692 primaryHistogramCollection.m_hMomentumEfficiency->GetXaxis()->SetRangeUser(1., +5.5);
693 primaryHistogramCollection.m_hMomentumEfficiency->GetXaxis()->SetTitle("True Momentum [GeV]");
694 primaryHistogramCollection.m_hMomentumEfficiency->GetYaxis()->SetRangeUser(0., +1.01);
695 primaryHistogramCollection.m_hMomentumEfficiency->GetYaxis()->SetTitle("Reconstruction Efficiency");
696 }
697
698 if (!primaryHistogramCollection.m_hCompleteness)
699 {
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");
704 }
705
706 if (!primaryHistogramCollection.m_hPurity)
707 {
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");
712 }
713
714 primaryHistogramCollection.m_hHitsAll->Fill(primaryResult.m_nMCHitsTotal);
715 primaryHistogramCollection.m_hMomentumAll->Fill(primaryResult.m_trueMomentum);
716
717 if ((primaryResult.m_nPfoMatches > 0) &&
718 (!parameters.m_correctTrackShowerId || primaryResult.m_isCorrectParticleId))
719 {
720 primaryHistogramCollection.m_hHitsEfficiency->Fill(primaryResult.m_nMCHitsTotal);
721 primaryHistogramCollection.m_hMomentumEfficiency->Fill(primaryResult.m_trueMomentum);
722 primaryHistogramCollection.m_hCompleteness->Fill(primaryResult.m_bestMatchCompleteness);
723 primaryHistogramCollection.m_hPurity->Fill(primaryResult.m_bestMatchPurity);
724 }
725}
726
727//------------------------------------------------------------------------------------------------------------------------------------------
728
729void ProcessHistogramCollections(const InteractionPrimaryHistogramMap &interactionPrimaryHistogramMap)
730{
731 for (InteractionPrimaryHistogramMap::const_iterator iter = interactionPrimaryHistogramMap.begin(), iterEnd = interactionPrimaryHistogramMap.end(); iter != iterEnd; ++iter)
732 {
733 const InteractionType interactionType(iter->first);
734 const PrimaryHistogramMap &primaryHistogramMap(iter->second);
735
736 for (PrimaryHistogramMap::const_iterator hIter = primaryHistogramMap.begin(), hIterEnd = primaryHistogramMap.end(); hIter != hIterEnd; ++hIter)
737 {
738 const ExpectedPrimary expectedPrimary(hIter->first);
739 const PrimaryHistogramCollection &primaryHistogramCollection(hIter->second);
740
741 for (int n = -1; n <= primaryHistogramCollection.m_hHitsEfficiency->GetXaxis()->GetNbins(); ++n)
742 {
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;
747 primaryHistogramCollection.m_hHitsEfficiency->SetBinContent(n + 1, efficiency);
748 primaryHistogramCollection.m_hHitsEfficiency->SetBinError(n + 1, error);
749 }
750
751 for (int n = -1; n <= primaryHistogramCollection.m_hMomentumEfficiency->GetXaxis()->GetNbins(); ++n)
752 {
753 const float found = primaryHistogramCollection.m_hMomentumEfficiency->GetBinContent(n + 1);
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;
757 primaryHistogramCollection.m_hMomentumEfficiency->SetBinContent(n + 1, efficiency);
758 primaryHistogramCollection.m_hMomentumEfficiency->SetBinError(n + 1, error);
759 }
760
761 primaryHistogramCollection.m_hCompleteness->Scale(1. / static_cast<double>(primaryHistogramCollection.m_hCompleteness->GetEntries()));
762 primaryHistogramCollection.m_hPurity->Scale(1. / static_cast<double>(primaryHistogramCollection.m_hPurity->GetEntries()));
763 }
764 }
765}
766
767//------------------------------------------------------------------------------------------------------------------------------------------
768
769std::string ToString(const ExpectedPrimary expectedPrimary)
770{
771 switch (expectedPrimary)
772 {
773 case MUON : return "MUON";
774 case ELECTRON : return "ELECTRON";
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";
785 case OTHER_PRIMARY: return "OTHER_PRIMARY";
786 default: return "UNKNOWN";
787 }
788}
789
790//------------------------------------------------------------------------------------------------------------------------------------------
791
792std::string ToString(const InteractionType interactionType)
793{
794 switch (interactionType)
795 {
796 case CCQEL_MU: return "CCQEL_MU";
797 case CCQEL_MU_P: return "CCQEL_MU_P";
798 case CCQEL_MU_P_P: return "CCQEL_MU_P_P";
799 case CCQEL_MU_P_P_P: return "CCQEL_MU_P_P_P";
800 case CCQEL_MU_P_P_P_P: return "CCQEL_MU_P_P_P_P";
801 case CCQEL_MU_P_P_P_P_P: return "CCQEL_MU_P_P_P_P_P";
802 case CCQEL_E: return "CCQEL_E";
803 case CCQEL_E_P: return "CCQEL_E_P";
804 case CCQEL_E_P_P: return "CCQEL_E_P_P";
805 case CCQEL_E_P_P_P: return "CCQEL_E_P_P_P";
806 case CCQEL_E_P_P_P_P: return "CCQEL_E_P_P_P_P";
807 case CCQEL_E_P_P_P_P_P: return "CCQEL_E_P_P_P_P_P";
808 case NCQEL_P: return "NCQEL_P";
809 case NCQEL_P_P: return "NCQEL_P_P";
810 case NCQEL_P_P_P: return "NCQEL_P_P_P";
811 case NCQEL_P_P_P_P: return "NCQEL_P_P_P_P";
812 case NCQEL_P_P_P_P_P: return "NCQEL_P_P_P_P_P";
813 case CCRES_MU: return "CCRES_MU";
814 case CCRES_MU_P: return "CCRES_MU_P";
815 case CCRES_MU_P_P: return "CCRES_MU_P_P";
816 case CCRES_MU_P_P_P: return "CCRES_MU_P_P_P";
817 case CCRES_MU_P_P_P_P: return "CCRES_MU_P_P_P_P";
818 case CCRES_MU_P_P_P_P_P: return "CCRES_MU_P_P_P_P_P";
819 case CCRES_MU_PIPLUS: return "CCRES_MU_PIPLUS";
820 case CCRES_MU_P_PIPLUS: return "CCRES_MU_P_PIPLUS";
821 case CCRES_MU_P_P_PIPLUS: return "CCRES_MU_P_P_PIPLUS";
822 case CCRES_MU_P_P_P_PIPLUS: return "CCRES_MU_P_P_P_PIPLUS";
823 case CCRES_MU_P_P_P_P_PIPLUS: return "CCRES_MU_P_P_P_P_PIPLUS";
824 case CCRES_MU_P_P_P_P_P_PIPLUS: return "CCRES_MU_P_P_P_P_P_PIPLUS";
825 case CCRES_MU_PHOTON: return "CCRES_MU_PHOTON";
826 case CCRES_MU_P_PHOTON: return "CCRES_MU_P_PHOTON";
827 case CCRES_MU_P_P_PHOTON: return "CCRES_MU_P_P_PHOTON";
828 case CCRES_MU_P_P_P_PHOTON: return "CCRES_MU_P_P_P_PHOTON";
829 case CCRES_MU_P_P_P_P_PHOTON: return "CCRES_MU_P_P_P_P_PHOTON";
830 case CCRES_MU_P_P_P_P_P_PHOTON: return "CCRES_MU_P_P_P_P_P_PHOTON";
831 case CCRES_MU_PIZERO: return "CCRES_MU_PIZERO";
832 case CCRES_MU_P_PIZERO: return "CCRES_MU_P_PIZERO";
833 case CCRES_MU_P_P_PIZERO: return "CCRES_MU_P_P_PIZERO";
834 case CCRES_MU_P_P_P_PIZERO: return "CCRES_MU_P_P_P_PIZERO";
835 case CCRES_MU_P_P_P_P_PIZERO: return "CCRES_MU_P_P_P_P_PIZERO";
836 case CCRES_MU_P_P_P_P_P_PIZERO: return "CCRES_MU_P_P_P_P_P_PIZERO";
837 case CCRES_E: return "CCRES_E";
838 case CCRES_E_P: return "CCRES_E_P";
839 case CCRES_E_P_P: return "CCRES_E_P_P";
840 case CCRES_E_P_P_P: return "CCRES_E_P_P_P";
841 case CCRES_E_P_P_P_P: return "CCRES_E_P_P_P_P";
842 case CCRES_E_P_P_P_P_P: return "CCRES_E_P_P_P_P_P";
843 case CCRES_E_PIPLUS: return "CCRES_E_PIPLUS";
844 case CCRES_E_P_PIPLUS: return "CCRES_E_P_PIPLUS";
845 case CCRES_E_P_P_PIPLUS: return "CCRES_E_P_P_PIPLUS";
846 case CCRES_E_P_P_P_PIPLUS: return "CCRES_E_P_P_P_PIPLUS";
847 case CCRES_E_P_P_P_P_PIPLUS: return "CCRES_E_P_P_P_P_PIPLUS";
848 case CCRES_E_P_P_P_P_P_PIPLUS: return "CCRES_E_P_P_P_P_P_PIPLUS";
849 case CCRES_E_PHOTON: return "CCRES_E_PHOTON";
850 case CCRES_E_P_PHOTON: return "CCRES_E_P_PHOTON";
851 case CCRES_E_P_P_PHOTON: return "CCRES_E_P_P_PHOTON";
852 case CCRES_E_P_P_P_PHOTON: return "CCRES_E_P_P_P_PHOTON";
853 case CCRES_E_P_P_P_P_PHOTON: return "CCRES_E_P_P_P_P_PHOTON";
854 case CCRES_E_P_P_P_P_P_PHOTON: return "CCRES_E_P_P_P_P_P_PHOTON";
855 case CCRES_E_PIZERO: return "CCRES_E_PIZERO";
856 case CCRES_E_P_PIZERO: return "CCRES_E_P_PIZERO";
857 case CCRES_E_P_P_PIZERO: return "CCRES_E_P_P_PIZERO";
858 case CCRES_E_P_P_P_PIZERO: return "CCRES_E_P_P_P_PIZERO";
859 case CCRES_E_P_P_P_P_PIZERO: return "CCRES_E_P_P_P_P_PIZERO";
860 case CCRES_E_P_P_P_P_P_PIZERO: return "CCRES_E_P_P_P_P_P_PIZERO";
861 case NCRES_P: return "NCRES_P";
862 case NCRES_P_P: return "NCRES_P_P";
863 case NCRES_P_P_P: return "NCRES_P_P_P";
864 case NCRES_P_P_P_P: return "NCRES_P_P_P_P";
865 case NCRES_P_P_P_P_P: return "NCRES_P_P_P_P_P";
866 case NCRES_PIPLUS: return "NCRES_PIPLUS";
867 case NCRES_P_PIPLUS: return "NCRES_P_PIPLUS";
868 case NCRES_P_P_PIPLUS: return "NCRES_P_P_PIPLUS";
869 case NCRES_P_P_P_PIPLUS: return "NCRES_P_P_P_PIPLUS";
870 case NCRES_P_P_P_P_PIPLUS: return "NCRES_P_P_P_P_PIPLUS";
871 case NCRES_P_P_P_P_P_PIPLUS: return "NCRES_P_P_P_P_P_PIPLUS";
872 case NCRES_PIMINUS: return "NCRES_PIMINUS";
873 case NCRES_P_PIMINUS: return "NCRES_P_PIMINUS";
874 case NCRES_P_P_PIMINUS: return "NCRES_P_P_PIMINUS";
875 case NCRES_P_P_P_PIMINUS: return "NCRES_P_P_P_PIMINUS";
876 case NCRES_P_P_P_P_PIMINUS: return "NCRES_P_P_P_P_PIMINUS";
877 case NCRES_P_P_P_P_P_PIMINUS: return "NCRES_P_P_P_P_P_PIMINUS";
878 case NCRES_PHOTON: return "NCRES_PHOTON";
879 case NCRES_P_PHOTON: return "NCRES_P_PHOTON";
880 case NCRES_P_P_PHOTON: return "NCRES_P_P_PHOTON";
881 case NCRES_P_P_P_PHOTON: return "NCRES_P_P_P_PHOTON";
882 case NCRES_P_P_P_P_PHOTON: return "NCRES_P_P_P_P_PHOTON";
883 case NCRES_P_P_P_P_P_PHOTON: return "NCRES_P_P_P_P_P_PHOTON";
884 case NCRES_PIZERO: return "NCRES_PIZERO";
885 case NCRES_P_PIZERO: return "NCRES_P_PIZERO";
886 case NCRES_P_P_PIZERO: return "NCRES_P_P_PIZERO";
887 case NCRES_P_P_P_PIZERO: return "NCRES_P_P_P_PIZERO";
888 case NCRES_P_P_P_P_PIZERO: return "NCRES_P_P_P_P_PIZERO";
889 case NCRES_P_P_P_P_P_PIZERO: return "NCRES_P_P_P_P_P_PIZERO";
890 case CCDIS_MU: return "CCDIS_MU";
891 case CCDIS_MU_P: return "CCDIS_MU_P";
892 case CCDIS_MU_P_P: return "CCDIS_MU_P_P";
893 case CCDIS_MU_P_P_P: return "CCDIS_MU_P_P_P";
894 case CCDIS_MU_P_P_P_P: return "CCDIS_MU_P_P_P_P";
895 case CCDIS_MU_P_P_P_P_P: return "CCDIS_MU_P_P_P_P_P";
896 case CCDIS_MU_PIPLUS: return "CCDIS_MU_PIPLUS";
897 case CCDIS_MU_P_PIPLUS: return "CCDIS_MU_P_PIPLUS";
898 case CCDIS_MU_P_P_PIPLUS: return "CCDIS_MU_P_P_PIPLUS";
899 case CCDIS_MU_P_P_P_PIPLUS: return "CCDIS_MU_P_P_P_PIPLUS";
900 case CCDIS_MU_P_P_P_P_PIPLUS: return "CCDIS_MU_P_P_P_P_PIPLUS";
901 case CCDIS_MU_P_P_P_P_P_PIPLUS: return "CCDIS_MU_P_P_P_P_P_PIPLUS";
902 case CCDIS_MU_PHOTON: return "CCDIS_MU_PHOTON";
903 case CCDIS_MU_P_PHOTON: return "CCDIS_MU_P_PHOTON";
904 case CCDIS_MU_P_P_PHOTON: return "CCDIS_MU_P_P_PHOTON";
905 case CCDIS_MU_P_P_P_PHOTON: return "CCDIS_MU_P_P_P_PHOTON";
906 case CCDIS_MU_P_P_P_P_PHOTON: return "CCDIS_MU_P_P_P_P_PHOTON";
907 case CCDIS_MU_P_P_P_P_P_PHOTON: return "CCDIS_MU_P_P_P_P_P_PHOTON";
908 case CCDIS_MU_PIZERO: return "CCDIS_MU_PIZERO";
909 case CCDIS_MU_P_PIZERO: return "CCDIS_MU_P_PIZERO";
910 case CCDIS_MU_P_P_PIZERO: return "CCDIS_MU_P_P_PIZERO";
911 case CCDIS_MU_P_P_P_PIZERO: return "CCDIS_MU_P_P_P_PIZERO";
912 case CCDIS_MU_P_P_P_P_PIZERO: return "CCDIS_MU_P_P_P_P_PIZERO";
913 case CCDIS_MU_P_P_P_P_P_PIZERO: return "CCDIS_MU_P_P_P_P_P_PIZERO";
914 case NCDIS_P: return "NCDIS_P";
915 case NCDIS_P_P: return "NCDIS_P_P";
916 case NCDIS_P_P_P: return "NCDIS_P_P_P";
917 case NCDIS_P_P_P_P: return "NCDIS_P_P_P_P";
918 case NCDIS_P_P_P_P_P: return "NCDIS_P_P_P_P_P";
919 case NCDIS_PIPLUS: return "NCDIS_PIPLUS";
920 case NCDIS_P_PIPLUS: return "NCDIS_P_PIPLUS";
921 case NCDIS_P_P_PIPLUS: return "NCDIS_P_P_PIPLUS";
922 case NCDIS_P_P_P_PIPLUS: return "NCDIS_P_P_P_PIPLUS";
923 case NCDIS_P_P_P_P_PIPLUS: return "NCDIS_P_P_P_P_PIPLUS";
924 case NCDIS_P_P_P_P_P_PIPLUS: return "NCDIS_P_P_P_P_P_PIPLUS";
925 case NCDIS_PIMINUS: return "NCDIS_PIMINUS";
926 case NCDIS_P_PIMINUS: return "NCDIS_P_PIMINUS";
927 case NCDIS_P_P_PIMINUS: return "NCDIS_P_P_PIMINUS";
928 case NCDIS_P_P_P_PIMINUS: return "NCDIS_P_P_P_PIMINUS";
929 case NCDIS_P_P_P_P_PIMINUS: return "NCDIS_P_P_P_P_PIMINUS";
930 case NCDIS_P_P_P_P_P_PIMINUS: return "NCDIS_P_P_P_P_P_PIMINUS";
931 case NCDIS_PHOTON: return "NCDIS_PHOTON";
932 case NCDIS_P_PHOTON: return "NCDIS_P_PHOTON";
933 case NCDIS_P_P_PHOTON: return "NCDIS_P_P_PHOTON";
934 case NCDIS_P_P_P_PHOTON: return "NCDIS_P_P_P_PHOTON";
935 case NCDIS_P_P_P_P_PHOTON: return "NCDIS_P_P_P_P_PHOTON";
936 case NCDIS_P_P_P_P_P_PHOTON: return "NCDIS_P_P_P_P_P_PHOTON";
937 case NCDIS_PIZERO: return "NCDIS_PIZERO";
938 case NCDIS_P_PIZERO: return "NCDIS_P_PIZERO";
939 case NCDIS_P_P_PIZERO: return "NCDIS_P_P_PIZERO";
940 case NCDIS_P_P_P_PIZERO: return "NCDIS_P_P_P_PIZERO";
941 case NCDIS_P_P_P_P_PIZERO: return "NCDIS_P_P_P_P_PIZERO";
942 case NCDIS_P_P_P_P_P_PIZERO: return "NCDIS_P_P_P_P_P_PIZERO";
943 case CCCOH: return "CCCOH";
944 case NCCOH: return "NCCOH";
945 case COSMIC_RAY_MU: return "COSMIC_RAY_MU";
946 case COSMIC_RAY_P: return "COSMIC_RAY_P";
947 case COSMIC_RAY_E: return "COSMIC_RAY_E";
948 case COSMIC_RAY_PHOTON: return "COSMIC_RAY_PHOTON";
949 case COSMIC_RAY_OTHER: return "COSMIC_RAY_OTHER";
950 case BEAM_PARTICLE_MU: return "BEAM_PARTICLE_MU";
951 case BEAM_PARTICLE_P: return "BEAM_PARTICLE_P";
952 case BEAM_PARTICLE_E: return "BEAM_PARTICLE_E";
953 case BEAM_PARTICLE_PHOTON: return "BEAM_PARTICLE_PHOTON";
954 case BEAM_PARTICLE_PI_PLUS: return "BEAM_PARTICLE_PI_PLUS";
955 case BEAM_PARTICLE_PI_MINUS: return "BEAM_PARTICLE_PI_MINUS";
956 case BEAM_PARTICLE_KAON_PLUS: return "BEAM_PARTICLE_KAON_PLUS";
957 case BEAM_PARTICLE_KAON_MINUS: return "BEAM_PARTICLE_KAON_MINUS";
958 case BEAM_PARTICLE_OTHER: return "BEAM_PARTICLE_OTHER";
959 case OTHER_INTERACTION: return "OTHER_INTERACTION";
960 case ALL_INTERACTIONS: return "ALL_INTERACTIONS";
961 default: return "UNKNOWN";
962 }
963}
bool PassSBNDFiducialCut(const SimpleMCTarget &simpleMCTarget)
Whether a simple mc event passes sbnd fiducial cut, applied to target vertices.
Definition Validation.C:415
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.
Definition Validation.C:469
void DisplayInteractionCountingMap(const InteractionCountingMap &interactionCountingMap, const Parameters &parameters)
Print details to screen for a provided interaction type to counting map.
Definition Validation.C:484
void CountPfoMatches(const SimpleMCEvent &simpleMCEvent, const Parameters &parameters, InteractionCountingMap &interactionCountingMap, InteractionTargetResultMap &interactionTargetResultMap)
CountPfoMatches Relies on fact that primary list is sorted by number of true good hits.
Definition Validation.C:311
int ReadNextEvent(TChain *const pTChain, const int iEntry, SimpleMCEvent &simpleMCEvent, const Parameters &parameters)
Read the next event from the chain.
Definition Validation.C:55
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...
Definition Validation.C:433
bool PassUbooneFiducialCut(const SimpleMCTarget &simpleMCTarget)
Whether a simple mc event passes uboone fiducial cut, applied to target vertices.
Definition Validation.C:398
bool PassFiducialCut(const SimpleMCTarget &simpleMCTarget, const Parameters &parameters)
Whether a simple mc event passes the relevant fiducial cut, applied to target vertices.
Definition Validation.C:382
void DisplaySimpleMCEventMatches(const SimpleMCEvent &simpleMCEvent, const Parameters &parameters)
Print matching details to screen for a simple mc event.
Definition Validation.C:212
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...
Definition Validation.C:613
void ProcessHistogramCollections(const InteractionPrimaryHistogramMap &interactionPrimaryHistogramMap)
Process histograms stored in the provided map e.g. calculating final efficiencies,...
Definition Validation.C:729
void FillPrimaryHistogramCollection(const std::string &histPrefix, const Parameters &parameters, const PrimaryResult &primaryResult, PrimaryHistogramCollection &primaryHistogramCollection)
Fill histograms in the provided histogram collection, using information in the provided primary resul...
Definition Validation.C:655
void AnalyseInteractionTargetResultMap(const InteractionTargetResultMap &interactionTargetResultMap, const Parameters &parameters)
Opportunity to fill histograms, perform post-processing of information collected in main loop over nt...
Definition Validation.C:533
void Validation(const std::string &inputFiles, const Parameters &parameters)
Validation - Main entry point for analysis.
Definition Validation.C:18
std::string ToString(const ExpectedPrimary expectedPrimary)
Get a string representation of an interaction type.
Definition Validation.C:769
std::map< InteractionType, PrimaryHistogramMap > InteractionPrimaryHistogramMap
Definition Validation.h:529
std::map< ExpectedPrimary, CountingDetails > CountingMap
Definition Validation.h:432
std::map< ExpectedPrimary, PrimaryHistogramCollection > PrimaryHistogramMap
Definition Validation.h:528
std::map< InteractionType, CountingMap > InteractionCountingMap
Definition Validation.h:433
std::map< ExpectedPrimary, PrimaryResult > PrimaryResultMap
Definition Validation.h:458
std::map< InteractionType, TargetResultList > InteractionTargetResultMap
Definition Validation.h:482
std::vector< TargetResult > TargetResultList
Definition Validation.h:481
std::vector< float > FloatVector
Definition Validation.h:14
InteractionType
InteractionType enum.
Definition Validation.h:234
@ NCRES_P_P_P_P_PIMINUS
Definition Validation.h:315
@ CCRES_E_P_P_P_P_P_PIZERO
Definition Validation.h:299
@ NCRES_P_PIPLUS
Definition Validation.h:306
@ CCRES_MU_P_P_P_PIPLUS
Definition Validation.h:261
@ CCRES_MU_P
Definition Validation.h:253
@ CCDIS_MU_P_P_P_P
Definition Validation.h:333
@ CCRES_E_PIPLUS
Definition Validation.h:282
@ CCDIS_MU_P_PIZERO
Definition Validation.h:348
@ COSMIC_RAY_P
Definition Validation.h:385
@ NCDIS_P_P
Definition Validation.h:354
@ NCQEL_P_P
Definition Validation.h:248
@ CCDIS_MU_PIZERO
Definition Validation.h:347
@ CCRES_MU_P_P_P
Definition Validation.h:255
@ NCDIS_PIZERO
Definition Validation.h:376
@ BEAM_PARTICLE_E
Definition Validation.h:391
@ CCDIS_MU_PHOTON
Definition Validation.h:341
@ NCRES_P_PIMINUS
Definition Validation.h:312
@ CCQEL_E_P_P_P
Definition Validation.h:244
@ CCRES_E_P_P
Definition Validation.h:278
@ CCRES_E_PIZERO
Definition Validation.h:294
@ BEAM_PARTICLE_KAON_MINUS
Definition Validation.h:396
@ NCDIS_P_P_P_P_PIZERO
Definition Validation.h:380
@ NCQEL_P
Definition Validation.h:247
@ CCQEL_MU_P
Definition Validation.h:236
@ CCRES_MU_P_PIPLUS
Definition Validation.h:259
@ CCQEL_MU_P_P_P
Definition Validation.h:238
@ NCQEL_P_P_P_P_P
Definition Validation.h:251
@ CCDIS_MU_P_P_P_PIPLUS
Definition Validation.h:338
@ CCRES_E_P_P_P_P_P_PIPLUS
Definition Validation.h:287
@ CCRES_MU_P_P_P_PHOTON
Definition Validation.h:267
@ CCRES_MU_P_P_P_P_PIPLUS
Definition Validation.h:262
@ CCDIS_MU_P_P_P_PIZERO
Definition Validation.h:350
@ NCDIS_PIPLUS
Definition Validation.h:358
@ CCDIS_MU
Definition Validation.h:329
@ NCRES_P_P_P_P_PIPLUS
Definition Validation.h:309
@ CCRES_E_P_P_PIZERO
Definition Validation.h:296
@ CCRES_MU_P_PHOTON
Definition Validation.h:265
@ NCDIS_P_P_PIPLUS
Definition Validation.h:360
@ CCDIS_MU_P_P_PIZERO
Definition Validation.h:349
@ CCDIS_MU_P_P_P_P_PIPLUS
Definition Validation.h:339
@ NCDIS_P_P_P_P
Definition Validation.h:356
@ CCDIS_MU_P_PIPLUS
Definition Validation.h:336
@ CCDIS_MU_P_P_PHOTON
Definition Validation.h:343
@ BEAM_PARTICLE_PI_MINUS
Definition Validation.h:394
@ CCRES_MU_P_P_P_PIZERO
Definition Validation.h:273
@ CCRES_E_P_PIPLUS
Definition Validation.h:283
@ NCDIS_P_PHOTON
Definition Validation.h:371
@ NCRES_P_P_P_P_P_PHOTON
Definition Validation.h:322
@ BEAM_PARTICLE_MU
Definition Validation.h:389
@ NCRES_P_P_P
Definition Validation.h:302
@ CCRES_E_P_P_P_PIZERO
Definition Validation.h:297
@ CCRES_E_P_PIZERO
Definition Validation.h:295
@ NCDIS_P_P_P_PIMINUS
Definition Validation.h:367
@ COSMIC_RAY_PHOTON
Definition Validation.h:387
@ BEAM_PARTICLE_KAON_PLUS
Definition Validation.h:395
@ CCRES_MU_P_P_P_P_P_PHOTON
Definition Validation.h:269
@ CCDIS_MU_P
Definition Validation.h:330
@ CCRES_MU_PHOTON
Definition Validation.h:264
@ BEAM_PARTICLE_P
Definition Validation.h:390
@ CCDIS_MU_P_P_P_P_P_PIPLUS
Definition Validation.h:340
@ NCDIS_PHOTON
Definition Validation.h:370
@ CCQEL_MU_P_P_P_P_P
Definition Validation.h:240
@ CCQEL_E_P_P_P_P_P
Definition Validation.h:246
@ NCDIS_P_P_P_P_PIPLUS
Definition Validation.h:362
@ NCDIS_P
Definition Validation.h:353
@ CCRES_MU_PIZERO
Definition Validation.h:270
@ NCRES_P
Definition Validation.h:300
@ NCDIS_P_P_P_P_P_PIZERO
Definition Validation.h:381
@ CCDIS_MU_P_P_P
Definition Validation.h:332
@ CCRES_MU_P_P
Definition Validation.h:254
@ CCDIS_MU_P_PHOTON
Definition Validation.h:342
@ NCDIS_P_P_P_P_P_PHOTON
Definition Validation.h:375
@ CCRES_MU_P_P_P_P_P_PIPLUS
Definition Validation.h:263
@ CCDIS_MU_P_P_P_P_PHOTON
Definition Validation.h:345
@ NCRES_P_P_PIMINUS
Definition Validation.h:313
@ CCRES_E_P_P_P_P_P_PHOTON
Definition Validation.h:293
@ COSMIC_RAY_E
Definition Validation.h:386
@ NCQEL_P_P_P_P
Definition Validation.h:250
@ CCRES_MU_PIPLUS
Definition Validation.h:258
@ CCRES_MU_P_P_PIZERO
Definition Validation.h:272
@ CCRES_E_P_P_PIPLUS
Definition Validation.h:284
@ NCDIS_P_P_P_P_PIMINUS
Definition Validation.h:368
@ NCRES_P_P_PIPLUS
Definition Validation.h:307
@ NCRES_P_P_P_PIMINUS
Definition Validation.h:314
@ CCRES_MU_P_P_P_P
Definition Validation.h:256
@ NCRES_PIMINUS
Definition Validation.h:311
@ NCRES_P_P_P_P
Definition Validation.h:303
@ NCRES_PIZERO
Definition Validation.h:323
@ NCRES_P_P_P_P_P
Definition Validation.h:304
@ NCDIS_P_P_P_P_P
Definition Validation.h:357
@ CCRES_MU_P_P_P_P_P_PIZERO
Definition Validation.h:275
@ CCRES_MU
Definition Validation.h:252
@ CCRES_E_P_P_P
Definition Validation.h:279
@ NCRES_P_P
Definition Validation.h:301
@ CCQEL_MU_P_P_P_P
Definition Validation.h:239
@ CCDIS_MU_P_P_P_PHOTON
Definition Validation.h:344
@ CCRES_E_P_P_P_P_PHOTON
Definition Validation.h:292
@ NCDIS_P_P_P_P_PHOTON
Definition Validation.h:374
@ NCRES_P_P_P_P_PHOTON
Definition Validation.h:321
@ NCRES_PIPLUS
Definition Validation.h:305
@ NCDIS_P_P_P_PHOTON
Definition Validation.h:373
@ CCRES_E_P_P_P_PIPLUS
Definition Validation.h:285
@ CCRES_MU_P_P_P_P_PHOTON
Definition Validation.h:268
@ NCDIS_P_PIPLUS
Definition Validation.h:359
@ NCCOH
Definition Validation.h:383
@ NCRES_P_P_P_P_P_PIPLUS
Definition Validation.h:310
@ NCRES_P_PHOTON
Definition Validation.h:318
@ NCRES_P_P_PHOTON
Definition Validation.h:319
@ CCDIS_MU_P_P
Definition Validation.h:331
@ CCQEL_MU
Definition Validation.h:235
@ NCDIS_P_PIMINUS
Definition Validation.h:365
@ NCDIS_P_P_P_P_P_PIMINUS
Definition Validation.h:369
@ CCDIS_MU_P_P_PIPLUS
Definition Validation.h:337
@ CCDIS_MU_P_P_P_P_P_PIZERO
Definition Validation.h:352
@ NCDIS_P_P_PHOTON
Definition Validation.h:372
@ NCRES_P_P_P_PHOTON
Definition Validation.h:320
@ NCRES_P_P_P_PIZERO
Definition Validation.h:326
@ CCRES_MU_P_PIZERO
Definition Validation.h:271
@ BEAM_PARTICLE_OTHER
Definition Validation.h:397
@ NCRES_P_P_P_P_P_PIZERO
Definition Validation.h:328
@ CCRES_E_P_P_P_P_PIPLUS
Definition Validation.h:286
@ CCRES_MU_P_P_PIPLUS
Definition Validation.h:260
@ CCQEL_E
Definition Validation.h:241
@ NCDIS_P_P_PIMINUS
Definition Validation.h:366
@ NCDIS_P_P_P_PIPLUS
Definition Validation.h:361
@ BEAM_PARTICLE_PHOTON
Definition Validation.h:392
@ NCRES_PHOTON
Definition Validation.h:317
@ CCRES_MU_P_P_P_P_P
Definition Validation.h:257
@ CCQEL_E_P
Definition Validation.h:242
@ NCDIS_P_P_P_PIZERO
Definition Validation.h:379
@ CCDIS_MU_P_P_P_P_P
Definition Validation.h:334
@ CCRES_E_P_PHOTON
Definition Validation.h:289
@ CCQEL_E_P_P_P_P
Definition Validation.h:245
@ CCRES_MU_P_P_P_P_PIZERO
Definition Validation.h:274
@ CCDIS_MU_P_P_P_P_PIZERO
Definition Validation.h:351
@ NCRES_P_P_P_P_P_PIMINUS
Definition Validation.h:316
@ NCRES_P_PIZERO
Definition Validation.h:324
@ CCQEL_MU_P_P
Definition Validation.h:237
@ NCRES_P_P_P_PIPLUS
Definition Validation.h:308
@ CCDIS_MU_PIPLUS
Definition Validation.h:335
@ CCDIS_MU_P_P_P_P_P_PHOTON
Definition Validation.h:346
@ NCQEL_P_P_P
Definition Validation.h:249
@ CCRES_E_P_P_P_P_PIZERO
Definition Validation.h:298
@ ALL_INTERACTIONS
Definition Validation.h:399
@ OTHER_INTERACTION
Definition Validation.h:398
@ BEAM_PARTICLE_PI_PLUS
Definition Validation.h:393
@ NCRES_P_P_P_P_PIZERO
Definition Validation.h:327
@ CCRES_E_P_P_PHOTON
Definition Validation.h:290
@ CCRES_E_P_P_P_PHOTON
Definition Validation.h:291
@ NCDIS_P_P_PIZERO
Definition Validation.h:378
@ CCQEL_E_P_P
Definition Validation.h:243
@ NCDIS_P_P_P_P_P_PIPLUS
Definition Validation.h:363
@ NCDIS_P_P_P
Definition Validation.h:355
@ CCRES_MU_P_P_PHOTON
Definition Validation.h:266
@ CCRES_E_P_P_P_P_P
Definition Validation.h:281
@ CCRES_E
Definition Validation.h:276
@ COSMIC_RAY_OTHER
Definition Validation.h:388
@ NCDIS_P_PIZERO
Definition Validation.h:377
@ CCCOH
Definition Validation.h:382
@ CCRES_E_P_P_P_P
Definition Validation.h:280
@ NCRES_P_P_PIZERO
Definition Validation.h:325
@ COSMIC_RAY_MU
Definition Validation.h:384
@ NCDIS_PIMINUS
Definition Validation.h:364
@ CCRES_E_PHOTON
Definition Validation.h:288
@ CCRES_E_P
Definition Validation.h:277
std::vector< int > IntVector
Definition Validation.h:13
ExpectedPrimary
ExpectedPrimary enum.
Definition Validation.h:203
@ NEUTRON
Definition Validation.h:213
@ PHOTON2
Definition Validation.h:215
@ PROTON4
Definition Validation.h:209
@ PROTON5
Definition Validation.h:210
@ PROTON3
Definition Validation.h:208
@ PHOTON1
Definition Validation.h:214
@ PROTON1
Definition Validation.h:206
@ ELECTRON
Definition Validation.h:205
@ PIMINUS
Definition Validation.h:212
@ PIPLUS
Definition Validation.h:211
@ PROTON2
Definition Validation.h:207
@ OTHER_PRIMARY
Definition Validation.h:216
@ MUON
Definition Validation.h:204
std::vector< SimpleMCPrimary > SimpleMCPrimaryList
Definition Validation.h:129
std::map< InteractionType, TargetHistogramCollection > InteractionTargetHistogramMap
Definition Validation.h:505
CountingDetails class.
Definition Validation.h:417
unsigned int m_nTotal
The total number of occurences.
Definition Validation.h:424
unsigned int m_correctId
The number of times the mc primary particle id was correct.
Definition Validation.h:429
unsigned int m_nMatch0
The number of times the mc primary has 0 pfo matches.
Definition Validation.h:425
unsigned int m_nMatch2
The number of times the mc primary has 2 pfo matches.
Definition Validation.h:427
unsigned int m_nMatch3Plus
The number of times the mc primary has 3 or more pfo matches.
Definition Validation.h:428
unsigned int m_nMatch1
The number of times the mc primary has 1 pfo matches.
Definition Validation.h:426
Parameters class.
Definition Validation.h:20
bool m_applyUbooneFiducialCut
Whether to apply uboone fiducial volume cut to true neutrino vertex position.
Definition Validation.h:30
bool m_triggeredBeamOnly
Whether to only consider triggered beam particles.
Definition Validation.h:36
std::string m_eventFileName
File name to which to write list of correct events.
Definition Validation.h:39
std::string m_histPrefix
Histogram name prefix.
Definition Validation.h:37
float m_vertexXCorrection
The vertex x correction, added to reported mc neutrino endpoint x value, in cm.
Definition Validation.h:33
std::string m_mapFileName
File name to which to write output ascii tables, etc.
Definition Validation.h:38
int m_nEventsToProcess
The number of events to process.
Definition Validation.h:29
bool m_displayMatchedEvents
Whether to display matching results for individual events.
Definition Validation.h:27
bool m_testBeamMode
Whether running in test beam mode.
Definition Validation.h:35
int m_skipEvents
The number of events to skip.
Definition Validation.h:28
bool m_histogramOutput
Whether to produce output histograms.
Definition Validation.h:34
bool m_correctTrackShowerId
Whether to demand that pfos are correctly flagged as tracks or showers.
Definition Validation.h:32
bool m_applySBNDFiducialCut
Whether to apply sbnd fiducial volume cut to true neutrino vertex position.
Definition Validation.h:31
PrimaryHistogramCollection class.
Definition Validation.h:513
TH1F * m_hMomentumAll
The number of primaries vs momentum histogram.
Definition Validation.h:522
TH1F * m_hHitsEfficiency
The primary efficiency vs number of hits histogram.
Definition Validation.h:521
TH1F * m_hPurity
The primary (best match) purity histogram.
Definition Validation.h:525
TH1F * m_hMomentumEfficiency
The primary efficiency vs momentum histogram.
Definition Validation.h:523
TH1F * m_hHitsAll
The number of primaries vs number of hits histogram.
Definition Validation.h:520
TH1F * m_hCompleteness
The primary (best match) completeness histogram.
Definition Validation.h:524
PrimaryResult class.
Definition Validation.h:441
bool m_isCorrectParticleId
Whether the best matched pfo has the correct particle id.
Definition Validation.h:454
unsigned int m_nPfoMatches
The total number of pfo matches for a given primary.
Definition Validation.h:448
unsigned int m_nMCHitsTotal
The number of hits in the mc primary.
Definition Validation.h:449
unsigned int m_nBestMatchRecoHitsTotal
The number of hits in the best matched pfo.
Definition Validation.h:451
float m_bestMatchPurity
The purity of the best matched pfo.
Definition Validation.h:453
float m_trueMomentum
The true momentum of the mc primary.
Definition Validation.h:455
float m_bestMatchCompleteness
The completeness of the best matched pfo.
Definition Validation.h:452
unsigned int m_nBestMatchSharedHitsTotal
The number of hits shared by the mc primary and the best matched pfo.
Definition Validation.h:450
SimpleMCEvent class.
Definition Validation.h:181
int m_eventNumber
The event number.
Definition Validation.h:189
int m_fileIdentifier
The file identifier.
Definition Validation.h:188
SimpleMCTargetList m_mcTargetList
The list of mc targets.
Definition Validation.h:192
int m_nMCTargets
The number of mc targets.
Definition Validation.h:191
SimpleMCPrimary class.
Definition Validation.h:93
SimpleThreeVector m_momentum
The momentum.
Definition Validation.h:103
int m_nMCHitsV
The number of v mc hits.
Definition Validation.h:108
int m_pdgCode
The pdg code.
Definition Validation.h:101
int m_bestMatchPfoIsTestBeam
Whether best match pfo is reconstructed as a test beam particle.
Definition Validation.h:118
int m_bestMatchPfoNHitsW
The best match pfo number of w pfo hits.
Definition Validation.h:122
int m_primaryId
The identifier.
Definition Validation.h:100
int m_bestMatchPfoNSharedHitsW
The best match pfo number of w matched hits.
Definition Validation.h:126
int m_bestMatchPfoIsRecoNu
Whether best match pfo is reconstructed as part of a neutrino hierarchy.
Definition Validation.h:116
int m_bestMatchPfoRecoNuId
The identifier of the associated reco neutrino (if part of a neutrino hierarchy)
Definition Validation.h:117
int m_nPrimaryMatchedPfos
The number of matched pfos.
Definition Validation.h:111
int m_bestMatchPfoNSharedHitsU
The best match pfo number of u matched hits.
Definition Validation.h:124
int m_bestMatchPfoNHitsTotal
The best match pfo total number of pfo hits.
Definition Validation.h:119
int m_nMCHitsTotal
The total number of mc hits.
Definition Validation.h:106
SimpleThreeVector m_vertex
The vertex.
Definition Validation.h:104
int m_bestMatchPfoNHitsU
The best match pfo number of u pfo hits.
Definition Validation.h:120
int m_nMCHitsU
The number of u mc hits.
Definition Validation.h:107
int m_bestMatchPfoNHitsV
The best match pfo number of v pfo hits.
Definition Validation.h:121
SimpleThreeVector m_endpoint
The endpoint.
Definition Validation.h:105
int m_bestMatchPfoId
The best match pfo identifier.
Definition Validation.h:114
int m_bestMatchPfoNSharedHitsV
The best match pfo number of v matched hits.
Definition Validation.h:125
int m_nPrimaryMatchedCRPfos
The number of matched cr pfos.
Definition Validation.h:113
float m_energy
The energy.
Definition Validation.h:102
int m_nMCHitsW
The number of w mc hits.
Definition Validation.h:109
int m_bestMatchPfoPdgCode
The best match pfo pdg code.
Definition Validation.h:115
int m_nPrimaryMatchedNuPfos
The number of matched nu pfos.
Definition Validation.h:112
int m_bestMatchPfoNSharedHitsTotal
The best match pfo total number of matched hits.
Definition Validation.h:123
SimpleMCTarget class.
Definition Validation.h:137
int m_isCosmicRay
Whether the target is a cosmic ray.
Definition Validation.h:148
int m_isFakeCR
Whether the target was reconstructed as a fake cosmic ray.
Definition Validation.h:157
int m_isNeutrino
Whether the target is a neutrino.
Definition Validation.h:146
int m_nTargetNuMatches
The number of neutrino pfo matches to the target.
Definition Validation.h:163
int m_isLost
Whether the target was lost (not reconstructed)
Definition Validation.h:160
int m_nTargetGoodNuMatches
The number of good neutrino pfo matches to the target (all from same parent neutrino)
Definition Validation.h:165
int m_isFakeNu
Whether the target was reconstructed as a fake neutrino.
Definition Validation.h:156
SimpleMCPrimaryList m_mcPrimaryList
The list of mc primaries.
Definition Validation.h:170
int m_nTargetNuLosses
The number of neutrino primaries with no matches.
Definition Validation.h:167
int m_isCorrectNu
Whether the target was correctly reconstructed as a neutrino.
Definition Validation.h:153
int m_interactionType
The target interaction type.
Definition Validation.h:144
int m_nTargetPrimaries
The number of target mc primaries.
Definition Validation.h:169
int m_isCorrectTB
Whether the target was correctly reconstructed as a beam particle.
Definition Validation.h:154
int m_nTargetMatches
The number of pfo matches to the target.
Definition Validation.h:162
int m_nTargetCRMatches
The number of cosmic ray pfo matches to the target.
Definition Validation.h:164
SimpleThreeVector m_targetVertex
The target vertex position.
Definition Validation.h:150
int m_nTargetNuSplits
The number of split neutrino pfo matches to the target (from different parent neutrinos)
Definition Validation.h:166
int m_isSplitNu
Whether the target was reconstructed as a split neutrino.
Definition Validation.h:158
int m_isCorrectCR
Whether the target was correctly reconstructed as a cosmic ray.
Definition Validation.h:155
SimpleThreeVector m_recoVertex
The reco vertex position, if available.
Definition Validation.h:151
int m_mcNuanceCode
The target nuance code.
Definition Validation.h:145
int m_isBeamParticle
Whether the target is a beam particle.
Definition Validation.h:147
int m_isSplitCR
Whether the target was reconstructed as a split cosmic ray.
Definition Validation.h:159
float m_x
The x value.
Definition Validation.h:64
float m_y
The y value.
Definition Validation.h:65
float m_z
The z value.
Definition Validation.h:66
TargetHistogramCollection class.
Definition Validation.h:492
TH1F * m_hVtxDeltaY
The vtx delta y histogram.
Definition Validation.h:500
TH1F * m_hVtxDeltaX
The vtx delta x histogram.
Definition Validation.h:499
TH1F * m_hVtxDeltaR
The vtx delta r histogram.
Definition Validation.h:502
TH1F * m_hVtxDeltaZ
The vtx delta z histogram.
Definition Validation.h:501
TargetResult class.
Definition Validation.h:466
bool m_hasRecoVertex
Whether a reco vertex is matched to the target.
Definition Validation.h:476
bool m_isCorrect
Whether the target is reconstructed correctly.
Definition Validation.h:475
PrimaryResultMap m_primaryResultMap
The primary result map.
Definition Validation.h:478
int m_fileIdentifier
The file identifier.
Definition Validation.h:473
int m_eventNumber
The event number.
Definition Validation.h:474
SimpleThreeVector m_vertexOffset
The offset between the reco and true target vertices.
Definition Validation.h:477