LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
GlobalAsymmetryFeatureTool.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
13 
14 using namespace pandora;
15 
16 namespace lar_content
17 {
18 
19 GlobalAsymmetryFeatureTool::GlobalAsymmetryFeatureTool() :
20  m_maxAsymmetryDistance(5.f)
21 {
22 }
23 
24 //------------------------------------------------------------------------------------------------------------------------------------------
25 
26 void GlobalAsymmetryFeatureTool::Run(LArMvaHelper::MvaFeatureVector &featureVector, const VertexSelectionBaseAlgorithm *const pAlgorithm, const Vertex * const pVertex,
29 {
30  if (PandoraContentApi::GetSettings(*pAlgorithm)->ShouldDisplayAlgorithmInfo())
31  std::cout << "----> Running Algorithm Tool: " << this->GetInstanceName() << ", " << this->GetType() << std::endl;
32 
33  float globalAsymmetry(0.f);
34 
35  globalAsymmetry += this->GetGlobalAsymmetryForView(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), TPC_VIEW_U),
36  slidingFitDataListMap.at(TPC_VIEW_U));
37 
38  globalAsymmetry += this->GetGlobalAsymmetryForView(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), TPC_VIEW_V),
39  slidingFitDataListMap.at(TPC_VIEW_V));
40 
41  globalAsymmetry += this->GetGlobalAsymmetryForView(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), TPC_VIEW_W),
42  slidingFitDataListMap.at(TPC_VIEW_W));
43 
44  featureVector.push_back(globalAsymmetry);
45 }
46 
47 //------------------------------------------------------------------------------------------------------------------------------------------
48 
49 float GlobalAsymmetryFeatureTool::GetGlobalAsymmetryForView(const CartesianVector &vertexPosition2D,
50  const VertexSelectionBaseAlgorithm::SlidingFitDataList &slidingFitDataList) const
51 {
52  bool useEnergy(true);
53  CartesianVector energyWeightedDirectionSum(0.f, 0.f, 0.f), hitWeightedDirectionSum(0.f, 0.f, 0.f);
54 
55  for (const VertexSelectionBaseAlgorithm::SlidingFitData &slidingFitData : slidingFitDataList)
56  {
57  const Cluster *const pCluster(slidingFitData.GetCluster());
58 
59  if (pCluster->GetElectromagneticEnergy() < std::numeric_limits<float>::epsilon())
60  useEnergy = false;
61 
62  const CartesianVector vertexToMinLayer(slidingFitData.GetMinLayerPosition() - vertexPosition2D);
63  const CartesianVector vertexToMaxLayer(slidingFitData.GetMaxLayerPosition() - vertexPosition2D);
64 
65  const bool minLayerClosest(vertexToMinLayer.GetMagnitudeSquared() < vertexToMaxLayer.GetMagnitudeSquared());
66  const CartesianVector &clusterDirection((minLayerClosest) ? slidingFitData.GetMinLayerDirection() : slidingFitData.GetMaxLayerDirection());
67 
68  if ((LArClusterHelper::GetClosestDistance(vertexPosition2D, pCluster) < m_maxAsymmetryDistance))
69  {
70  this->IncrementAsymmetryParameters(pCluster->GetElectromagneticEnergy(), clusterDirection, energyWeightedDirectionSum);
71  this->IncrementAsymmetryParameters(static_cast<float>(pCluster->GetNCaloHits()), clusterDirection, hitWeightedDirectionSum);
72  }
73  }
74 
75  const CartesianVector &localWeightedDirectionSum(useEnergy ? energyWeightedDirectionSum : hitWeightedDirectionSum);
76 
77  if (localWeightedDirectionSum.GetMagnitudeSquared() < std::numeric_limits<float>::epsilon())
78  return 0.f;
79 
80  return this->CalculateGlobalAsymmetry(useEnergy, vertexPosition2D, slidingFitDataList, localWeightedDirectionSum);
81 }
82 
83 //------------------------------------------------------------------------------------------------------------------------------------------
84 
85 void GlobalAsymmetryFeatureTool::IncrementAsymmetryParameters(const float weight, const CartesianVector &clusterDirection,
86  CartesianVector &localWeightedDirectionSum) const
87 {
88  // If the new axis direction is at an angle of greater than 90 deg to the current axis direction, flip it 180 degs.
89  CartesianVector newDirection(clusterDirection);
90 
91  if (localWeightedDirectionSum.GetMagnitudeSquared() > std::numeric_limits<float>::epsilon())
92  {
93  if (localWeightedDirectionSum.GetCosOpeningAngle(clusterDirection) < 0.f)
94  newDirection *= -1.f;
95  }
96 
97  localWeightedDirectionSum += newDirection * weight;
98 }
99 
100 //------------------------------------------------------------------------------------------------------------------------------------------
101 
102 float GlobalAsymmetryFeatureTool::CalculateGlobalAsymmetry(const bool useEnergyMetrics, const CartesianVector &vertexPosition2D,
103  const VertexSelectionBaseAlgorithm::SlidingFitDataList &slidingFitDataList, const CartesianVector &localWeightedDirectionSum) const
104 {
105  // Project every hit onto local event axis direction and record side of the projected vtx position on which it falls
106  float beforeVtxHitEnergy(0.f), afterVtxHitEnergy(0.f);
107  unsigned int beforeVtxHitCount(0), afterVtxHitCount(0);
108 
109  const CartesianVector localWeightedDirection(localWeightedDirectionSum.GetUnitVector());
110  const float evtProjectedVtxPos(vertexPosition2D.GetDotProduct(localWeightedDirection));
111 
112  for (const VertexSelectionBaseAlgorithm::SlidingFitData &slidingFitData : slidingFitDataList)
113  {
114  const Cluster * const pCluster(slidingFitData.GetCluster());
115 
116  CaloHitList caloHitList;
117  pCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
118 
119  CaloHitVector caloHitVector(caloHitList.begin(), caloHitList.end());
120  std::sort(caloHitVector.begin(), caloHitVector.end(), LArClusterHelper::SortHitsByPosition);
121 
122  for (const CaloHit *const pCaloHit : caloHitVector)
123  {
124  if (pCaloHit->GetPositionVector().GetDotProduct(localWeightedDirection) < evtProjectedVtxPos)
125  {
126  beforeVtxHitEnergy += pCaloHit->GetElectromagneticEnergy();
127  ++beforeVtxHitCount;
128  }
129 
130  else
131  {
132  afterVtxHitEnergy += pCaloHit->GetElectromagneticEnergy();
133  ++afterVtxHitCount;
134  }
135  }
136  }
137 
138  // Use energy metrics if possible, otherwise fall back on hit counting.
139  const float totHitEnergy(beforeVtxHitEnergy + afterVtxHitEnergy);
140  const unsigned int totHitCount(beforeVtxHitCount + afterVtxHitCount);
141 
142  if (useEnergyMetrics)
143  return std::fabs((afterVtxHitEnergy - beforeVtxHitEnergy)) / totHitEnergy;
144 
145  if (0 == totHitCount)
146  throw StatusCodeException(STATUS_CODE_FAILURE);
147 
148  return std::fabs((static_cast<float>(afterVtxHitCount) - static_cast<float>(beforeVtxHitCount))) / static_cast<float>(totHitCount);
149 }
150 
151 //------------------------------------------------------------------------------------------------------------------------------------------
152 
153 StatusCode GlobalAsymmetryFeatureTool::ReadSettings(const TiXmlHandle xmlHandle)
154 {
155  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
156  "MaxAsymmetryDistance", m_maxAsymmetryDistance));
157 
158  return STATUS_CODE_SUCCESS;
159 }
160 
161 } // namespace lar_content
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
MvaTypes::MvaFeatureVector MvaFeatureVector
Definition: LArMvaHelper.h:58
static pandora::CartesianVector ProjectPosition(const pandora::Pandora &pandora, const pandora::CartesianVector &position3D, const pandora::HitType view)
Project 3D position into a given 2D view.
void IncrementAsymmetryParameters(const float weight, const pandora::CartesianVector &clusterDirection, pandora::CartesianVector &localWeightedDirectionSum) const
Increment the asymmetry parameters.
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
std::map< pandora::HitType, const ShowerClusterList > ShowerClusterListMap
Map of shower cluster lists for passing to tools.
static bool SortHitsByPosition(const pandora::CaloHit *const pLhs, const pandora::CaloHit *const pRhs)
Sort calo hits by their position (use Z, followed by X, followed by Y)
Header file for the cluster helper class.
std::map< pandora::HitType, const pandora::ClusterList & > ClusterListMap
Map array of cluster lists for passing to tools.
float CalculateGlobalAsymmetry(const bool useEnergyMetrics, const pandora::CartesianVector &vertexPosition2D, const VertexSelectionBaseAlgorithm::SlidingFitDataList &slidingFitDataList, const pandora::CartesianVector &localWeightedDirectionSum) const
Calculate the global asymmetry feature.
void Run(LArMvaHelper::MvaFeatureVector &featureVector, const VertexSelectionBaseAlgorithm *const pAlgorithm, const pandora::Vertex *const pVertex, const VertexSelectionBaseAlgorithm::SlidingFitDataListMap &slidingFitDataListMap, const VertexSelectionBaseAlgorithm::ClusterListMap &, const VertexSelectionBaseAlgorithm::KDTreeMap &, const VertexSelectionBaseAlgorithm::ShowerClusterListMap &, const float, float &)
Run the tool.
double weight
Definition: plottest35.C:25
float GetGlobalAsymmetryForView(const pandora::CartesianVector &vertexPosition2D, const VertexSelectionBaseAlgorithm::SlidingFitDataList &slidingFitDataList) const
Get the global asymmetry feature for a given view.
std::map< pandora::HitType, const SlidingFitDataList > SlidingFitDataListMap
Map of sliding fit data lists for passing to tools.
std::map< pandora::HitType, const std::reference_wrapper< HitKDTree2D > > KDTreeMap
Map array of hit kd trees for passing to tools.
float m_maxAsymmetryDistance
The max distance between cluster (any hit) and vertex to calculate asymmetry score.
Header file for the global asymmetry feature tool class.
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.