LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArMonitoringHelper.cc
Go to the documentation of this file.
1 
9 #include "Pandora/PdgTable.h"
10 
11 #include "Objects/CaloHit.h"
12 #include "Objects/MCParticle.h"
13 #include "Objects/ParticleFlowObject.h"
14 
15 #include "Helpers/MCParticleHelper.h"
16 
20 
21 #include <algorithm>
22 
23 using namespace pandora;
24 
25 namespace lar_content
26 {
27 
28 unsigned int LArMonitoringHelper::CountHitsByType(const HitType hitType, const CaloHitList &caloHitList)
29 {
30  unsigned int nHitsOfSpecifiedType(0);
31 
32  for (const CaloHit *const pCaloHit : caloHitList)
33  {
34  if (hitType == pCaloHit->GetHitType())
35  ++nHitsOfSpecifiedType;
36  }
37 
38  return nHitsOfSpecifiedType;
39 }
40 
41 //------------------------------------------------------------------------------------------------------------------------------------------
42 
43 void LArMonitoringHelper::GetOrderedMCParticleVector(
44  const LArMCParticleHelper::MCContributionMapVector &selectedMCParticleToGoodHitsMaps, MCParticleVector &orderedMCParticleVector)
45 {
46  for (const LArMCParticleHelper::MCContributionMap &mcParticleToGoodHitsMap : selectedMCParticleToGoodHitsMaps)
47  {
48  if (mcParticleToGoodHitsMap.empty())
49  continue;
50 
51  // Copy map contents to vector it can be sorted
52  std::vector<LArMCParticleHelper::MCParticleCaloHitListPair> mcParticleToGoodHitsVect;
53  std::copy(mcParticleToGoodHitsMap.begin(), mcParticleToGoodHitsMap.end(), std::back_inserter(mcParticleToGoodHitsVect));
54 
55  // Sort by number of hits descending
56  std::sort(mcParticleToGoodHitsVect.begin(), mcParticleToGoodHitsVect.end(),
58  {
59  // Neutrinos, then beam particles, then cosmic rays
60  const bool isANuFinalState(LArMCParticleHelper::IsBeamNeutrinoFinalState(a.first)),
61  isBNuFinalState(LArMCParticleHelper::IsBeamNeutrinoFinalState(b.first));
62 
63  if (isANuFinalState != isBNuFinalState)
64  return isANuFinalState;
65 
66  const bool isABeamParticle(LArMCParticleHelper::IsBeamParticle(a.first)),
67  isBBeamParticle(LArMCParticleHelper::IsBeamParticle(b.first));
68 
69  if (isABeamParticle != isBBeamParticle)
70  return isABeamParticle;
71 
72  // Then sort by numbers of true hits
73  if (a.second.size() != b.second.size())
74  return (a.second.size() > b.second.size());
75 
76  // Default to normal MCParticle sorting
77  return LArMCParticleHelper::SortByMomentum(a.first, b.first);
78  });
79 
80  for (const LArMCParticleHelper::MCParticleCaloHitListPair &mcParticleCaloHitPair : mcParticleToGoodHitsVect)
81  orderedMCParticleVector.push_back(mcParticleCaloHitPair.first);
82  }
83 
84  // Check that all elements of the vector are unique
85  const unsigned int nMCParticles(orderedMCParticleVector.size());
86  if (std::distance(orderedMCParticleVector.begin(), std::unique(orderedMCParticleVector.begin(), orderedMCParticleVector.end())) != nMCParticles)
87  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
88 }
89 
90 //------------------------------------------------------------------------------------------------------------------------------------------
91 
92 void LArMonitoringHelper::GetOrderedPfoVector(const LArMCParticleHelper::PfoContributionMap &pfoToReconstructable2DHitsMap, pandora::PfoVector &orderedPfoVector)
93 {
94  // Copy map contents to vector it can be sorted
95  std::vector<LArMCParticleHelper::PfoCaloHitListPair> pfoToReconstructable2DHitsVect;
96  std::copy(pfoToReconstructable2DHitsMap.begin(), pfoToReconstructable2DHitsMap.end(), std::back_inserter(pfoToReconstructable2DHitsVect));
97 
98  // Sort by number of hits descending putting neutrino final states first
99  std::sort(pfoToReconstructable2DHitsVect.begin(), pfoToReconstructable2DHitsVect.end(),
101  {
102  // Neutrinos before cosmic rays
103  const bool isANuFinalState(LArPfoHelper::IsNeutrinoFinalState(a.first)), isBNuFinalState(LArPfoHelper::IsNeutrinoFinalState(b.first));
104 
105  if (isANuFinalState != isBNuFinalState)
106  return isANuFinalState;
107 
108  if (a.second.size() != b.second.size())
109  return (a.second.size() > b.second.size());
110 
111  // Default to normal pfo sorting
112  return LArPfoHelper::SortByNHits(a.first, b.first);
113  });
114 
115  for (const LArMCParticleHelper::PfoCaloHitListPair &pfoCaloHitPair : pfoToReconstructable2DHitsVect)
116  orderedPfoVector.push_back(pfoCaloHitPair.first);
117 
118  // Check that all elements of the vector are unique
119  const unsigned int nPfos(orderedPfoVector.size());
120  if (std::distance(orderedPfoVector.begin(), std::unique(orderedPfoVector.begin(), orderedPfoVector.end())) != nPfos)
121  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
122 }
123 
124 //------------------------------------------------------------------------------------------------------------------------------------------
125 
126 void LArMonitoringHelper::PrintMCParticleTable(
127  const LArMCParticleHelper::MCContributionMap &selectedMCParticleToGoodHitsMap, const MCParticleVector &orderedMCParticleVector)
128 {
129  if (selectedMCParticleToGoodHitsMap.empty())
130  {
131  std::cout << "No MCParticles supplied." << std::endl;
132  return;
133  }
134 
135  LArFormattingHelper::Table table({"ID", "NUANCE", "TYPE", "", "E", "dist", "", "nGoodHits", "U", "V", "W"});
136 
137  unsigned int usedParticleCount(0);
138  for (unsigned int id = 0; id < orderedMCParticleVector.size(); ++id)
139  {
140  const MCParticle *const pMCParticle(orderedMCParticleVector.at(id));
141 
142  LArMCParticleHelper::MCContributionMap::const_iterator it = selectedMCParticleToGoodHitsMap.find(pMCParticle);
143  if (selectedMCParticleToGoodHitsMap.end() == it)
144  continue; // ATTN MCParticles in selectedMCParticleToGoodHitsMap may be a subset of orderedMCParticleVector
145 
146  // ATTN enumerate from 1 to match event validation algorithm
147  table.AddElement(id + 1);
148  table.AddElement(LArMCParticleHelper::GetNuanceCode(pMCParticle));
149  table.AddElement(PdgTable::GetParticleName(pMCParticle->GetParticleId()));
150 
151  table.AddElement(pMCParticle->GetEnergy());
152  table.AddElement((pMCParticle->GetEndpoint() - pMCParticle->GetVertex()).GetMagnitude());
153 
154  table.AddElement(it->second.size());
155  table.AddElement(LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, it->second));
156  table.AddElement(LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, it->second));
157  table.AddElement(LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, it->second));
158 
159  usedParticleCount++;
160  }
161 
162  // Check every MCParticle in selectedMCParticleToGoodHitsMap has been printed
163  if (usedParticleCount != selectedMCParticleToGoodHitsMap.size())
164  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
165 
166  table.Print();
167 }
168 
169 //------------------------------------------------------------------------------------------------------------------------------------------
170 
171 void LArMonitoringHelper::PrintPfoTable(const LArMCParticleHelper::PfoContributionMap &pfoToReconstructable2DHitsMap, const PfoVector &orderedPfoVector)
172 {
173  if (pfoToReconstructable2DHitsMap.empty())
174  {
175  std::cout << "No Pfos supplied." << std::endl;
176  return;
177  }
178 
179  LArFormattingHelper::Table table({"ID", "PID", "Is Nu FS", "", "nGoodHits", "U", "V", "W"});
180 
181  for (unsigned int id = 0; id < orderedPfoVector.size(); ++id)
182  {
183  const ParticleFlowObject *const pPfo(orderedPfoVector.at(id));
184 
185  LArMCParticleHelper::PfoContributionMap::const_iterator it = pfoToReconstructable2DHitsMap.find(pPfo);
186  if (pfoToReconstructable2DHitsMap.end() == it)
187  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
188 
189  // ATTN enumerate from 1 to match event validation algorithm
190  table.AddElement(id + 1);
191  table.AddElement(pPfo->GetParticleId());
192  table.AddElement(LArPfoHelper::IsNeutrinoFinalState(pPfo));
193 
194  table.AddElement(it->second.size());
195  table.AddElement(LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, it->second));
196  table.AddElement(LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, it->second));
197  table.AddElement(LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, it->second));
198  }
199 
200  table.Print();
201 }
202 
203 //------------------------------------------------------------------------------------------------------------------------------------------
204 
205 void LArMonitoringHelper::PrintMatchingTable(const PfoVector &orderedPfoVector, const MCParticleVector &orderedMCParticleVector,
206  const LArMCParticleHelper::MCParticleToPfoHitSharingMap &mcParticleToPfoHitSharingMap, const unsigned int nMatches)
207 {
208  if (orderedPfoVector.empty())
209  {
210  std::cout << "No Pfos supplied." << std::endl;
211  return;
212  }
213 
214  if (orderedMCParticleVector.empty())
215  {
216  std::cout << "No MCParticles supplied." << std::endl;
217  return;
218  }
219 
220  // Get the maximum number of MCParticle to Pfos matches that need to be shown
221  unsigned int maxMatches(0);
222  for (const auto &entry : mcParticleToPfoHitSharingMap)
223  maxMatches = std::max(static_cast<unsigned int>(entry.second.size()), maxMatches);
224 
225  const bool showOthersColumn(maxMatches > nMatches);
226  const unsigned int nMatchesToShow(std::min(maxMatches, nMatches));
227 
228  // Set up the table headers
229  std::vector<std::string> tableHeaders({"MCParticle", ""});
230  for (unsigned int i = 0; i < nMatchesToShow; ++i)
231  {
232  tableHeaders.push_back("");
233  tableHeaders.push_back("Pfo");
234  tableHeaders.push_back("nSharedHits");
235  }
236 
237  if (showOthersColumn)
238  {
239  tableHeaders.push_back("");
240  tableHeaders.push_back("");
241  tableHeaders.push_back("nOtherPfos");
242  tableHeaders.push_back("nSharedHits");
243  }
244 
245  LArFormattingHelper::Table table(tableHeaders);
246 
247  // Make a new row for each MCParticle
248  for (unsigned int mcParticleId = 0; mcParticleId < orderedMCParticleVector.size(); ++mcParticleId)
249  {
250  const MCParticle *const pMCParticle(orderedMCParticleVector.at(mcParticleId));
251  LArMCParticleHelper::MCParticleToPfoHitSharingMap::const_iterator it = mcParticleToPfoHitSharingMap.find(pMCParticle);
252  LArMCParticleHelper::PfoToSharedHitsVector pfoToSharedHitsVector;
253 
254  if (it != mcParticleToPfoHitSharingMap.end())
255  pfoToSharedHitsVector = it->second;
256 
257  const LArFormattingHelper::Color mcCol(LArMCParticleHelper::IsBeamNeutrinoFinalState(pMCParticle)
258  ? LArFormattingHelper::LIGHT_GREEN
259  : (LArMCParticleHelper::IsBeamParticle(pMCParticle) ? LArFormattingHelper::LIGHT_BLUE : LArFormattingHelper::LIGHT_RED));
260 
261  // ATTN enumerate from 1 to match event validation algorithm
262  table.AddElement(mcParticleId + 1, LArFormattingHelper::REGULAR, mcCol);
263 
264  // Get the matched Pfos
265  unsigned int nPfosShown(0);
266  unsigned int nOtherHits(0);
267  for (const auto &pfoNSharedHitsPair : pfoToSharedHitsVector)
268  {
269  for (unsigned int pfoId = 0; pfoId < orderedPfoVector.size(); ++pfoId)
270  {
271  if (pfoNSharedHitsPair.first != orderedPfoVector.at(pfoId))
272  continue;
273 
274  if (nPfosShown < nMatchesToShow)
275  {
276  // ATTN enumerate from 1 to match event validation algorithm
277  const LArFormattingHelper::Color pfoCol(
278  LArPfoHelper::IsNeutrinoFinalState(pfoNSharedHitsPair.first) ? LArFormattingHelper::LIGHT_GREEN : LArFormattingHelper::LIGHT_RED);
279  table.AddElement(pfoId + 1, LArFormattingHelper::REGULAR, pfoCol);
280  table.AddElement(pfoNSharedHitsPair.second.size(), LArFormattingHelper::REGULAR, pfoCol);
281  nPfosShown++;
282  }
283  else
284  {
285  nOtherHits += pfoNSharedHitsPair.second.size();
286  }
287  break;
288  }
289  }
290 
291  // Pad the rest of the row with empty entries
292  for (unsigned int i = 0; i < nMatchesToShow - nPfosShown; ++i)
293  {
294  table.AddElement("");
295  table.AddElement("");
296  }
297 
298  // Print any remaining matches
299  if (!showOthersColumn)
300  continue;
301 
302  if (nOtherHits != 0)
303  {
304  table.AddElement(pfoToSharedHitsVector.size() - nPfosShown);
305  table.AddElement(nOtherHits);
306  }
307  else
308  {
309  table.AddElement("");
310  table.AddElement("");
311  }
312  }
313  table.Print();
314 }
315 
316 } // namespace lar_content
Header file for the pfo helper class.
std::unordered_map< const pandora::MCParticle *, pandora::CaloHitList > MCContributionMap
std::pair< const pandora::MCParticle *, pandora::CaloHitList > MCParticleCaloHitListPair
std::map< const pandora::MCParticle *, PfoToSharedHitsVector > MCParticleToPfoHitSharingMap
std::vector< PfoCaloHitListPair > PfoToSharedHitsVector
intermediate_table::const_iterator const_iterator
Header file for the lar monitoring helper helper class.
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
void AddElement(const T &value, const Style style=REGULAR, const Color color=DEFAULT)
Add an element to the table into the next (non-separator) column.
std::vector< MCContributionMap > MCContributionMapVector
std::pair< const pandora::ParticleFlowObject *, pandora::CaloHitList > PfoCaloHitListPair
HitType
Definition: HitType.h:12
void Print() const
Print the table.
std::unordered_map< const pandora::ParticleFlowObject *, pandora::CaloHitList > PfoContributionMap