LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LayerSplittingAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
13 using namespace pandora;
14 
15 namespace lar_content
16 {
17 
18 LayerSplittingAlgorithm::LayerSplittingAlgorithm() :
19  m_minClusterLayers(20),
20  m_layerWindow(10),
21  m_maxScatterRms(0.35f),
22  m_maxScatterCosTheta(0.5f),
23  m_maxSlidingCosTheta(0.866f)
24 {
25 }
26 
27 //------------------------------------------------------------------------------------------------------------------------------------------
28 
29 StatusCode LayerSplittingAlgorithm::DivideCaloHits(const Cluster *const pCluster, CaloHitList &firstHitList, CaloHitList &secondHitList) const
30 {
31  unsigned int splitLayer(0);
32 
33  if (STATUS_CODE_SUCCESS == this->FindBestSplitLayer(pCluster, splitLayer))
34  return this->DivideCaloHits(pCluster, splitLayer, firstHitList, secondHitList);
35 
36  return STATUS_CODE_NOT_FOUND;
37 }
38 
39 //------------------------------------------------------------------------------------------------------------------------------------------
40 
41 StatusCode LayerSplittingAlgorithm::FindBestSplitLayer(const Cluster *const pCluster, unsigned int &splitLayer) const
42 {
43  const OrderedCaloHitList &orderedCaloHitList(pCluster->GetOrderedCaloHitList());
44 
45  if (orderedCaloHitList.size() < m_minClusterLayers)
46  return STATUS_CODE_NOT_FOUND;
47 
48  bool foundSplit(false);
49 
50  float bestCosTheta(1.f);
51  CartesianVector bestPosition(0.f, 0.f, 0.f);
52 
53  for (unsigned int iLayer = pCluster->GetInnerPseudoLayer() + 4; iLayer + 4 <= pCluster->GetOuterPseudoLayer(); ++iLayer)
54  {
55  if (orderedCaloHitList.find(iLayer) == orderedCaloHitList.end())
56  continue;
57 
58  unsigned int innerLayer((pCluster->GetInnerPseudoLayer() + m_layerWindow > iLayer) ? pCluster->GetInnerPseudoLayer() : iLayer - m_layerWindow);
59  unsigned int outerLayer((iLayer + m_layerWindow > pCluster->GetOuterPseudoLayer()) ? pCluster->GetOuterPseudoLayer() : iLayer + m_layerWindow);
60 
61  for (; innerLayer >= pCluster->GetInnerPseudoLayer(); --innerLayer)
62  {
63  if (orderedCaloHitList.find(innerLayer) != orderedCaloHitList.end())
64  break;
65  }
66 
67  for (; outerLayer <= pCluster->GetOuterPseudoLayer(); ++outerLayer)
68  {
69  if (orderedCaloHitList.find(outerLayer) != orderedCaloHitList.end())
70  break;
71  }
72 
73  const CartesianVector splitPosition(pCluster->GetCentroid(iLayer));
74  const CartesianVector innerPosition(pCluster->GetCentroid(innerLayer));
75  const CartesianVector outerPosition(pCluster->GetCentroid(outerLayer));
76 
77  const CartesianVector r1(innerPosition - splitPosition);
78  const CartesianVector r2(outerPosition - splitPosition);
79  const CartesianVector p1(r1.GetUnitVector());
80  const CartesianVector p2(r2.GetUnitVector());
81 
82  const float cosTheta(-p1.GetDotProduct(p2));
83  const float rms1(this->CalculateRms(pCluster, innerLayer, iLayer));
84  const float rms2(this->CalculateRms(pCluster, outerLayer, iLayer));
85  const float rms(std::max(rms1, rms2));
86 
87  float rmsCut(std::numeric_limits<float>::max());
88 
89  if (cosTheta > 0.f)
90  {
91  rmsCut = m_maxScatterRms;
92 
93  if (cosTheta > m_maxScatterCosTheta)
94  {
95  rmsCut *= ((m_maxSlidingCosTheta > cosTheta) ? (m_maxSlidingCosTheta - cosTheta) / (m_maxSlidingCosTheta - m_maxScatterCosTheta) : 0.f);
96  }
97  }
98 
99  if (rms < rmsCut && cosTheta < bestCosTheta)
100  {
101  bestCosTheta = cosTheta;
102  bestPosition = splitPosition;
103 
104  splitLayer = iLayer;
105  foundSplit = true;
106  }
107  }
108 
109  if (!foundSplit)
110  return STATUS_CODE_NOT_FOUND;
111 
112  return STATUS_CODE_SUCCESS;
113 }
114 
115 //------------------------------------------------------------------------------------------------------------------------------------------
116 
117 float LayerSplittingAlgorithm::CalculateRms(const Cluster *const pCluster, const unsigned int &firstLayer, const unsigned int &secondLayer) const
118 {
119  const OrderedCaloHitList &orderedCaloHitList(pCluster->GetOrderedCaloHitList());
120 
121  const unsigned int innerLayer(std::min(firstLayer, secondLayer));
122  const unsigned int outerLayer(std::max(firstLayer, secondLayer));
123 
124  const CartesianVector innerPosition(pCluster->GetCentroid(innerLayer));
125  const CartesianVector outerPosition(pCluster->GetCentroid(outerLayer));
126  const CartesianVector predictedDirection((outerPosition - innerPosition).GetUnitVector());
127 
128  float totalChi2(0.f);
129  float totalLayers(0.f);
130 
131  for (unsigned int iLayer = innerLayer + 1; iLayer + 1 < outerLayer; ++iLayer)
132  {
133  if (orderedCaloHitList.find(iLayer) == orderedCaloHitList.end())
134  continue;
135 
136  const CartesianVector hitPosition(pCluster->GetCentroid(iLayer));
137  const CartesianVector predictedPosition(innerPosition + predictedDirection * predictedDirection.GetDotProduct(hitPosition - innerPosition));
138 
139  totalChi2 += (predictedPosition - hitPosition).GetMagnitudeSquared();
140  totalLayers += 1.f;
141  }
142 
143  if (totalLayers > 0.f)
144  return std::sqrt(totalChi2 / totalLayers);
145 
146  return 0.f;
147 }
148 
149 //------------------------------------------------------------------------------------------------------------------------------------------
150 
152  const Cluster *const pCluster, const unsigned int &splitLayer, CaloHitList &firstHitList, CaloHitList &secondHitList) const
153 {
154  const OrderedCaloHitList &orderedCaloHitList(pCluster->GetOrderedCaloHitList());
155 
156  for (OrderedCaloHitList::const_iterator iter = orderedCaloHitList.begin(); iter != orderedCaloHitList.end(); ++iter)
157  {
158  const unsigned int thisLayer(iter->first);
159 
160  for (CaloHitList::const_iterator hitIter = iter->second->begin(), hitIterEnd = iter->second->end(); hitIter != hitIterEnd; ++hitIter)
161  {
162  const CaloHit *const pCaloHit = *hitIter;
163 
164  if (thisLayer < splitLayer)
165  {
166  firstHitList.push_back(pCaloHit);
167  }
168  else
169  {
170  secondHitList.push_back(pCaloHit);
171  }
172  }
173  }
174 
175  if (firstHitList.empty() || secondHitList.empty())
176  return STATUS_CODE_NOT_FOUND;
177 
178  return STATUS_CODE_SUCCESS;
179 }
180 
181 //------------------------------------------------------------------------------------------------------------------------------------------
182 
183 StatusCode LayerSplittingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
184 {
185  PANDORA_RETURN_RESULT_IF_AND_IF(
186  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterLayers", m_minClusterLayers));
187 
188  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "LayerWindow", m_layerWindow));
189 
190  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxScatterRms", m_maxScatterRms));
191 
192  PANDORA_RETURN_RESULT_IF_AND_IF(
193  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxScatterCosTheta", m_maxScatterCosTheta));
194 
195  PANDORA_RETURN_RESULT_IF_AND_IF(
196  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxSlidingCosTheta", m_maxSlidingCosTheta));
197 
198  return ClusterSplittingAlgorithm::ReadSettings(xmlHandle);
199 }
200 
201 } // namespace lar_content
pandora::StatusCode FindBestSplitLayer(const pandora::Cluster *const pCluster, unsigned int &splitLayer) const
Find the best layer for splitting the cluster.
intermediate_table::const_iterator const_iterator
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
TFile f
Definition: plotHisto.C:6
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
float CalculateRms(const pandora::Cluster *const pCluster, const unsigned int &firstLayer, const unsigned int &secondLayer) const
Calculate rms deviation of cluster centroids between two extremal layers.
pandora::StatusCode DivideCaloHits(const pandora::Cluster *const pCluster, pandora::CaloHitList &firstCaloHitList, pandora::CaloHitList &secondCaloHitList) const
Divide calo hits in a cluster into two lists, each associated with a separate fragment cluster...