LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PreProcessingAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
14 
16 
17 using namespace pandora;
18 
19 namespace lar_content
20 {
21 
22 PreProcessingAlgorithm::PreProcessingAlgorithm() :
23  m_mipEquivalentCut(std::numeric_limits<float>::epsilon()),
24  m_minCellLengthScale(std::numeric_limits<float>::epsilon()),
25  m_maxCellLengthScale(3.f),
26  m_searchRegion1D(0.1f),
27  m_maxEventHits(std::numeric_limits<unsigned int>::max()),
28  m_onlyAvailableCaloHits(true),
29  m_inputCaloHitListName("Input")
30 {
31 }
32 
33 //------------------------------------------------------------------------------------------------------------------------------------------
34 
36 {
37  m_processedHits.clear();
38  return STATUS_CODE_SUCCESS;
39 }
40 
41 //------------------------------------------------------------------------------------------------------------------------------------------
42 
44 {
45  if (!this->GetPandora().GetSettings()->SingleHitTypeClusteringMode())
46  {
47  std::cout << "PreProcessingAlgorithm: expect Pandora to be configured in SingleHitTypeClusteringMode." << std::endl;
48  return STATUS_CODE_FAILURE;
49  }
50 
51  try
52  {
53  this->ProcessCaloHits();
54  }
55  catch (const StatusCodeException &statusCodeException)
56  {
57  if (STATUS_CODE_OUT_OF_RANGE == statusCodeException.GetStatusCode())
58  {
59  std::cout << "PreProcessingAlgorithm: Excessive number of hits in event, skipping the reconstruction" << std::endl;
61  }
62  }
63 
65  {
66  if (STATUS_CODE_SUCCESS != PandoraContentApi::ReplaceCurrentList<CaloHit>(*this, m_currentCaloHitListReplacement))
67  {
68  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
69  std::cout << "PreProcessingAlgorithm: could not replace current calo hit list with list named: " << m_currentCaloHitListReplacement
70  << std::endl;
71  }
72  }
73 
74  return STATUS_CODE_SUCCESS;
75 }
76 
77 //------------------------------------------------------------------------------------------------------------------------------------------
78 
80 {
81  const CaloHitList *pCaloHitList(nullptr);
82  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*this, m_inputCaloHitListName, pCaloHitList));
83 
84  if (pCaloHitList->empty())
85  return;
86 
87  if (pCaloHitList->size() > m_maxEventHits)
88  throw StatusCodeException(STATUS_CODE_OUT_OF_RANGE);
89 
90  CaloHitList selectedCaloHitListU, selectedCaloHitListV, selectedCaloHitListW;
91 
92  for (const CaloHit *const pCaloHit : *pCaloHitList)
93  {
94  if (m_processedHits.count(pCaloHit))
95  continue;
96 
97  (void)m_processedHits.insert(pCaloHit);
98 
99  if (m_onlyAvailableCaloHits && !PandoraContentApi::IsAvailable(*this, pCaloHit))
100  continue;
101 
102  if (pCaloHit->GetMipEquivalentEnergy() < m_mipEquivalentCut)
103  continue;
104 
105  if (pCaloHit->GetInputEnergy() < std::numeric_limits<float>::epsilon())
106  {
107  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
108  std::cout << "PreProcessingAlgorithm: found a hit with zero energy, will remove it" << std::endl;
109 
110  continue;
111  }
112 
113  if ((pCaloHit->GetCellLengthScale() < m_minCellLengthScale) || (pCaloHit->GetCellLengthScale() > m_maxCellLengthScale))
114  {
115  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
116  {
117  std::cout << "PreProcessingAlgorithm: found a hit with extent " << pCaloHit->GetCellLengthScale() << ", require ("
118  << m_minCellLengthScale << " - " << m_maxCellLengthScale << "), will remove it" << std::endl;
119  }
120 
121  continue;
122  }
123 
124  if (TPC_VIEW_U == pCaloHit->GetHitType())
125  {
126  selectedCaloHitListU.push_back(pCaloHit);
127  }
128  else if (TPC_VIEW_V == pCaloHit->GetHitType())
129  {
130  selectedCaloHitListV.push_back(pCaloHit);
131  }
132  else if (TPC_VIEW_W == pCaloHit->GetHitType())
133  {
134  selectedCaloHitListW.push_back(pCaloHit);
135  }
136  }
137 
138  CaloHitList filteredCaloHitListU, filteredCaloHitListV, filteredCaloHitListW;
139  this->GetFilteredCaloHitList(selectedCaloHitListU, filteredCaloHitListU);
140  this->GetFilteredCaloHitList(selectedCaloHitListV, filteredCaloHitListV);
141  this->GetFilteredCaloHitList(selectedCaloHitListW, filteredCaloHitListW);
142 
143  CaloHitList filteredInputList;
144  filteredInputList.insert(filteredInputList.end(), filteredCaloHitListU.begin(), filteredCaloHitListU.end());
145  filteredInputList.insert(filteredInputList.end(), filteredCaloHitListV.begin(), filteredCaloHitListV.end());
146  filteredInputList.insert(filteredInputList.end(), filteredCaloHitListW.begin(), filteredCaloHitListW.end());
147 
148  if (!filteredInputList.empty() && !m_filteredCaloHitListName.empty())
149  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*this, filteredInputList, m_filteredCaloHitListName));
150 
151  if (!filteredCaloHitListU.empty() && !m_outputCaloHitListNameU.empty())
152  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*this, filteredCaloHitListU, m_outputCaloHitListNameU));
153 
154  if (!filteredCaloHitListV.empty() && !m_outputCaloHitListNameV.empty())
155  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*this, filteredCaloHitListV, m_outputCaloHitListNameV));
156 
157  if (!filteredCaloHitListW.empty() && !m_outputCaloHitListNameW.empty())
158  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*this, filteredCaloHitListW, m_outputCaloHitListNameW));
159 }
160 
161 //------------------------------------------------------------------------------------------------------------------------------------------
162 
164 {
165  CaloHitList emptyList;
166 
167  try
168  {
169  if (!m_filteredCaloHitListName.empty())
170  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*this, emptyList, m_filteredCaloHitListName));
171 
172  if (!m_outputCaloHitListNameU.empty())
173  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*this, emptyList, m_outputCaloHitListNameU));
174 
175  if (!m_outputCaloHitListNameV.empty())
176  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*this, emptyList, m_outputCaloHitListNameV));
177 
178  if (!m_outputCaloHitListNameW.empty())
179  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*this, emptyList, m_outputCaloHitListNameW));
180  }
181  catch (...)
182  {
183  std::cout << "PreProcessingAlgorithm: Unable to create void calo hit lists" << std::endl;
184  }
185 }
186 
187 //------------------------------------------------------------------------------------------------------------------------------------------
188 
189 void PreProcessingAlgorithm::GetFilteredCaloHitList(const CaloHitList &inputList, CaloHitList &outputList)
190 {
191  HitKDTree2D kdTree;
192  HitKDNode2DList hitKDNode2DList;
193 
194  KDTreeBox hitsBoundingRegion2D = fill_and_bound_2d_kd_tree(inputList, hitKDNode2DList);
195  kdTree.build(hitKDNode2DList, hitsBoundingRegion2D);
196 
197  // Remove hits that are in the same physical location!
198  for (const CaloHit *const pCaloHit1 : inputList)
199  {
200  bool isUnique(true);
202 
203  HitKDNode2DList found;
204  kdTree.search(searchRegionHits, found);
205 
206  for (const auto &hit : found)
207  {
208  const CaloHit *const pCaloHit2(hit.data);
209 
210  if (pCaloHit1 == pCaloHit2)
211  continue;
212 
213  const float displacementSquared((pCaloHit2->GetPositionVector() - pCaloHit1->GetPositionVector()).GetMagnitudeSquared());
214 
215  if (displacementSquared < std::numeric_limits<float>::epsilon())
216  {
217  const float deltaMip(pCaloHit2->GetMipEquivalentEnergy() > pCaloHit1->GetMipEquivalentEnergy());
218 
219  if ((deltaMip > std::numeric_limits<float>::epsilon()) ||
220  ((std::fabs(deltaMip) < std::numeric_limits<float>::epsilon()) &&
221  (outputList.end() != std::find(outputList.begin(), outputList.end(), pCaloHit2))))
222  {
223  isUnique = false;
224  break;
225  }
226  }
227  }
228 
229  if (isUnique)
230  {
231  outputList.push_back(pCaloHit1);
232  }
233  else
234  {
235  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
236  std::cout << "PreProcessingAlgorithm: found two hits in same location, will remove lowest pulse height" << std::endl;
237  }
238  }
239 }
240 
241 //------------------------------------------------------------------------------------------------------------------------------------------
242 
243 StatusCode PreProcessingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
244 {
245  PANDORA_RETURN_RESULT_IF_AND_IF(
246  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MipEquivalentCut", m_mipEquivalentCut));
247 
248  PANDORA_RETURN_RESULT_IF_AND_IF(
249  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCellLengthScale", m_minCellLengthScale));
250 
251  PANDORA_RETURN_RESULT_IF_AND_IF(
252  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxCellLengthScale", m_maxCellLengthScale));
253 
254  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SearchRegion1D", m_searchRegion1D));
255 
256  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxEventHits", m_maxEventHits));
257 
258  PANDORA_RETURN_RESULT_IF_AND_IF(
259  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "OnlyAvailableCaloHits", m_onlyAvailableCaloHits));
260 
261  PANDORA_RETURN_RESULT_IF_AND_IF(
262  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "InputCaloHitListName", m_inputCaloHitListName));
263 
264  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputCaloHitListNameU", m_outputCaloHitListNameU));
265 
266  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputCaloHitListNameV", m_outputCaloHitListNameV));
267 
268  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputCaloHitListNameW", m_outputCaloHitListNameW));
269 
270  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "FilteredCaloHitListName", m_filteredCaloHitListName));
271 
272  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "CurrentCaloHitListReplacement", m_currentCaloHitListReplacement));
273 
274  return STATUS_CODE_SUCCESS;
275 }
276 
277 } // namespace lar_content
std::string m_outputCaloHitListNameV
The output calo hit list name for TPC_VIEW_V hits.
Header file for the kd tree linker algo template class.
void GetFilteredCaloHitList(const pandora::CaloHitList &inputList, pandora::CaloHitList &outputList)
Clean up the input CaloHitList.
Box structure used to define 2D field. It&#39;s used in KDTree building step to divide the detector space...
float m_maxCellLengthScale
The maximum length scale for calo hit.
float m_mipEquivalentCut
Minimum mip equivalent energy for calo hit.
STL namespace.
unsigned int m_maxEventHits
The maximum number of hits in an event to proceed with the reconstruction.
std::string m_currentCaloHitListReplacement
The name of the calo hit list to replace the current list (optional)
TFile f
Definition: plotHisto.C:6
std::string m_outputCaloHitListNameW
The output calo hit list name for TPC_VIEW_W hits.
Header file for the cluster helper class.
void build(std::vector< KDTreeNodeInfoT< DATA, DIM >> &eltList, const KDTreeBoxT< DIM > &region)
Build the KD tree from the "eltList" in the space define by "region".
std::vector< HitKDNode2D > HitKDNode2DList
std::string m_filteredCaloHitListName
The output calo hit list name for all U, V and W hits.
float m_minCellLengthScale
The minimum length scale for calo hit.
std::string m_outputCaloHitListNameU
The output calo hit list name for TPC_VIEW_U hits.
Detector simulation of raw signals on wires.
std::string m_inputCaloHitListName
The input calo hit list name.
void ProcessCaloHits()
Build separate CaloHitLists for each view.
KDTreeBox fill_and_bound_2d_kd_tree(const MANAGED_CONTAINER< const T * > &points, std::vector< KDTreeNodeInfoT< const T *, 2 >> &nodes)
fill_and_bound_2d_kd_tree
void PopulateVoidCaloHitLists() noexcept
Build empty calo hit lists.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
float m_searchRegion1D
Search region, applied to each dimension, for look-up from kd-trees.
Header file for the pre processing algorithm class.
pandora::CaloHitSet m_processedHits
The set of all previously processed calo hits.
KDTreeBox build_2d_kd_search_region(const pandora::CaloHit *const point, const float x_span, const float z_span)
build_2d_kd_search_region
void search(const KDTreeBoxT< DIM > &searchBox, std::vector< KDTreeNodeInfoT< DATA, DIM >> &resRecHitList)
Search in the KDTree for all points that would be contained in the given searchbox The founded points...
bool m_onlyAvailableCaloHits
Whether to only include available calo hits.