LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArStitchingHelper.cc
Go to the documentation of this file.
1 
9 #include "Managers/GeometryManager.h"
10 
12 
14 
15 #include <cmath>
16 #include <limits>
17 
18 using namespace pandora;
19 
20 namespace lar_content
21 {
22 
23 const LArTPC &LArStitchingHelper::FindClosestTPC(const Pandora &pandora, const LArTPC &inputTPC, const bool checkPositive)
24 {
25  if (!MultiPandoraApi::GetPandoraInstanceMap().count(&pandora))
26  {
27  std::cout << "LArStitchingHelper::FindClosestTPC - functionality only available to primary/master Pandora instance " << std::endl;
28  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
29  }
30 
31  const LArTPCMap &larTPCMap(pandora.GetGeometry()->GetLArTPCMap());
32 
33  const LArTPC *pClosestTPC(nullptr);
34  float closestSeparation(std::numeric_limits<float>::max());
35  const float maxDisplacement(30.f); // TODO: 30cm should be fine, but can we do better than a hard-coded number here?
36 
37  for (const LArTPCMap::value_type &mapEntry : larTPCMap)
38  {
39  const LArTPC &checkTPC(*(mapEntry.second));
40 
41  if (&inputTPC == &checkTPC)
42  continue;
43 
44  if (checkPositive != (checkTPC.GetCenterX() > inputTPC.GetCenterX()))
45  continue;
46 
47  const float deltaX(std::fabs(checkTPC.GetCenterX() - inputTPC.GetCenterX()));
48  const float deltaY(std::fabs(checkTPC.GetCenterY() - inputTPC.GetCenterY()));
49  const float deltaZ(std::fabs(checkTPC.GetCenterZ() - inputTPC.GetCenterZ()));
50 
51  if (deltaY > maxDisplacement || deltaZ > maxDisplacement)
52  continue;
53 
54  if (deltaX < closestSeparation)
55  {
56  closestSeparation = deltaX;
57  pClosestTPC = &checkTPC;
58  }
59  }
60 
61  if (!pClosestTPC)
62  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
63 
64  return (*pClosestTPC);
65 }
66 
67 //------------------------------------------------------------------------------------------------------------------------------------------
68 
69 bool LArStitchingHelper::CanTPCsBeStitched(const LArTPC &firstTPC, const LArTPC &secondTPC)
70 {
71  if (&firstTPC == &secondTPC)
72  return false;
73 
74  // ATTN: We assume that Pfos are crossing either an anode-anode boundary or a cathode-cathode boundary
75  if (firstTPC.IsDriftInPositiveX() == secondTPC.IsDriftInPositiveX())
76  return false;
77 
78  return LArStitchingHelper::AreTPCsAdjacent(firstTPC, secondTPC);
79 }
80 
81 //------------------------------------------------------------------------------------------------------------------------------------------
82 
83 bool LArStitchingHelper::AreTPCsAdjacent(const LArTPC &firstTPC, const LArTPC &secondTPC)
84 {
85  // Check the relative positions of the centres of each drift volume
86  const float maxDisplacement(30.f); // TODO: 30cm should be fine, but can we do better than a hard-coded number here?
87  const float widthX(0.5f * (firstTPC.GetWidthX() + secondTPC.GetWidthX()));
88  const float deltaX(std::fabs(firstTPC.GetCenterX() - secondTPC.GetCenterX()));
89  const float deltaY(std::fabs(firstTPC.GetCenterY() - secondTPC.GetCenterY()));
90  const float deltaZ(std::fabs(firstTPC.GetCenterZ() - secondTPC.GetCenterZ()));
91 
92  if (std::fabs(deltaX - widthX) > maxDisplacement || deltaY > maxDisplacement || deltaZ > maxDisplacement)
93  return false;
94 
95  return true;
96 }
97 
98 //------------------------------------------------------------------------------------------------------------------------------------------
99 
100 bool LArStitchingHelper::AreTPCsAdjacent(const Pandora &pandora, const LArTPC &firstTPC, const LArTPC &secondTPC)
101 {
102  // Check if first tpc is just upstream of second tpc
103  try
104  {
105  const LArTPC &firstTPCCheck(LArStitchingHelper::FindClosestTPC(pandora, secondTPC, true));
106  const LArTPC &secondTPCCheck(LArStitchingHelper::FindClosestTPC(pandora, firstTPC, false));
107 
108  if ((&firstTPCCheck == &firstTPC) && (&secondTPCCheck == &secondTPC))
109  return true;
110  }
111  catch (pandora::StatusCodeException &)
112  {
113  }
114 
115  // Check if second tpc is just upstream of first tpc
116  try
117  {
118  const LArTPC &firstTPCCheck(LArStitchingHelper::FindClosestTPC(pandora, secondTPC, false));
119  const LArTPC &secondTPCCheck(LArStitchingHelper::FindClosestTPC(pandora, firstTPC, true));
120 
121  if ((&firstTPCCheck == &firstTPC) && (&secondTPCCheck == &secondTPC))
122  return true;
123  }
124  catch (pandora::StatusCodeException &)
125  {
126  }
127 
128  return false;
129 }
130 
131 //------------------------------------------------------------------------------------------------------------------------------------------
132 
133 float LArStitchingHelper::GetTPCBoundaryCenterX(const LArTPC &firstTPC, const LArTPC &secondTPC)
134 {
135  if (!LArStitchingHelper::AreTPCsAdjacent(firstTPC, secondTPC))
136  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
137 
138  if (firstTPC.GetCenterX() < secondTPC.GetCenterX())
139  {
140  return 0.5 * ((firstTPC.GetCenterX() + 0.5 * firstTPC.GetWidthX()) + (secondTPC.GetCenterX() - 0.5 * secondTPC.GetWidthX()));
141  }
142  else
143  {
144  return 0.5 * ((firstTPC.GetCenterX() - 0.5 * firstTPC.GetWidthX()) + (secondTPC.GetCenterX() + 0.5 * secondTPC.GetWidthX()));
145  }
146 }
147 
148 //------------------------------------------------------------------------------------------------------------------------------------------
149 
150 float LArStitchingHelper::GetTPCBoundaryWidthX(const LArTPC &firstTPC, const LArTPC &secondTPC)
151 {
152  if (!LArStitchingHelper::AreTPCsAdjacent(firstTPC, secondTPC))
153  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
154 
155  if (firstTPC.GetCenterX() < secondTPC.GetCenterX())
156  {
157  return ((secondTPC.GetCenterX() - 0.5 * secondTPC.GetWidthX()) - (firstTPC.GetCenterX() + 0.5 * firstTPC.GetWidthX()));
158  }
159  else
160  {
161  return ((firstTPC.GetCenterX() - 0.5 * firstTPC.GetWidthX()) - (secondTPC.GetCenterX() + 0.5 * secondTPC.GetWidthX()));
162  }
163 }
164 
165 //------------------------------------------------------------------------------------------------------------------------------------------
166 
167 float LArStitchingHelper::GetTPCDisplacement(const LArTPC &firstTPC, const LArTPC &secondTPC)
168 {
169  const float deltaX(firstTPC.GetCenterX() - secondTPC.GetCenterX());
170  const float deltaY(firstTPC.GetCenterY() - secondTPC.GetCenterY());
171  const float deltaZ(firstTPC.GetCenterZ() - secondTPC.GetCenterZ());
172 
173  return std::sqrt(deltaX * deltaX + deltaY * deltaY + deltaZ * deltaZ);
174 }
175 
176 //------------------------------------------------------------------------------------------------------------------------------------------
177 
178 void LArStitchingHelper::GetClosestVertices(const LArTPC &larTPC1, const LArTPC &larTPC2, const LArPointingCluster &pointingCluster1,
179  const LArPointingCluster &pointingCluster2, LArPointingCluster::Vertex &closestVertex1, LArPointingCluster::Vertex &closestVertex2)
180 {
181  if (&larTPC1 == &larTPC2)
182  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
183 
184  // Find the closest vertices based on relative X positions in tpc
185  const float dxVolume(larTPC2.GetCenterX() - larTPC1.GetCenterX());
186  const float dx1(pointingCluster1.GetOuterVertex().GetPosition().GetX() - pointingCluster1.GetInnerVertex().GetPosition().GetX());
187  const float dx2(pointingCluster2.GetOuterVertex().GetPosition().GetX() - pointingCluster2.GetInnerVertex().GetPosition().GetX());
188 
189  if (std::fabs(dx1) < std::numeric_limits<float>::epsilon() || std::fabs(dx2) < std::numeric_limits<float>::epsilon())
190  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
191 
192  const bool useInner1((dxVolume > 0.f) == (dx1 < 0.f)); // xVol1 - OUTER - INNER - | - xVol2 [xVol2-xVol1>0; xOuter1-xInner1<0]
193  const bool useInner2((dxVolume > 0.f) == (dx2 > 0.f)); // xVol1 - | - INNER - OUTER - xVol2 [xVol2-xVol1>0; xOuter2-xInner2>0]
194 
195  // Confirm that these really are the closest pair of vertices by checking the other possible pairs
196  const LArPointingCluster::Vertex &nearVertex1(useInner1 ? pointingCluster1.GetInnerVertex() : pointingCluster1.GetOuterVertex());
197  const LArPointingCluster::Vertex &nearVertex2(useInner2 ? pointingCluster2.GetInnerVertex() : pointingCluster2.GetOuterVertex());
198 
199  const LArPointingCluster::Vertex &farVertex1(useInner1 ? pointingCluster1.GetOuterVertex() : pointingCluster1.GetInnerVertex());
200  const LArPointingCluster::Vertex &farVertex2(useInner2 ? pointingCluster2.GetOuterVertex() : pointingCluster2.GetInnerVertex());
201 
202  const float dxNearNear(0.f);
203  const float dyNearNear(nearVertex1.GetPosition().GetY() - nearVertex2.GetPosition().GetY());
204  const float dzNearNear(nearVertex1.GetPosition().GetZ() - nearVertex2.GetPosition().GetZ());
205  const float drNearNearSquared(dxNearNear * dxNearNear + dyNearNear * dyNearNear + dzNearNear * dzNearNear);
206 
207  const float dxNearFar(std::fabs(dx2));
208  const float dyNearFar(nearVertex1.GetPosition().GetY() - farVertex2.GetPosition().GetY());
209  const float dzNearFar(nearVertex1.GetPosition().GetZ() - farVertex2.GetPosition().GetZ());
210  const float drNearFarSquared(dxNearFar * dxNearFar + dyNearFar * dyNearFar + dzNearFar * dzNearFar);
211 
212  const float dxFarNear(std::fabs(dx1));
213  const float dyFarNear(farVertex1.GetPosition().GetY() - nearVertex2.GetPosition().GetY());
214  const float dzFarNear(farVertex1.GetPosition().GetZ() - nearVertex2.GetPosition().GetZ());
215  const float drFarNearSquared(dxFarNear * dxFarNear + dyFarNear * dyFarNear + dzFarNear * dzFarNear);
216 
217  const float dxFarFar(std::fabs(dx1) + std::fabs(dx2));
218  const float dyFarFar(farVertex1.GetPosition().GetY() - farVertex2.GetPosition().GetY());
219  const float dzFarFar(farVertex1.GetPosition().GetZ() - farVertex2.GetPosition().GetZ());
220  const float drFarFarSquared(dxFarFar * dxFarFar + dyFarFar * dyFarFar + dzFarFar * dzFarFar);
221 
222  if (drNearNearSquared > std::min(drFarFarSquared, std::min(drNearFarSquared, drFarNearSquared)))
223  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
224 
225  closestVertex1 = nearVertex1;
226  closestVertex2 = nearVertex2;
227 }
228 
229 //------------------------------------------------------------------------------------------------------------------------------------------
230 
231 float LArStitchingHelper::CalculateX0(const LArTPC &firstTPC, const LArTPC &secondTPC, const LArPointingCluster::Vertex &firstVertex,
232  const LArPointingCluster::Vertex &secondVertex)
233 {
234  if (&firstTPC == &secondTPC)
235  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
236 
237  // ATTN: Assume that Pfos are crossing either an anode-anode boundary or a cathode-cathode boundary
238  if (firstTPC.IsDriftInPositiveX() == secondTPC.IsDriftInPositiveX())
239  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
240 
241  // Assume that Pfos have opposite direction components in x, and have some direction component in the y-z plane
242  const CartesianVector firstDirectionX(firstVertex.GetDirection().GetX(), 0.f, 0.f);
243  const CartesianVector secondDirectionX(secondVertex.GetDirection().GetX(), 0.f, 0.f);
244 
245  if (std::fabs(firstDirectionX.GetX()) > 1.f - std::numeric_limits<float>::epsilon() ||
246  std::fabs(secondDirectionX.GetX()) > 1.f - std::numeric_limits<float>::epsilon())
247  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
248 
249  if (-firstDirectionX.GetDotProduct(secondDirectionX) < std::numeric_limits<float>::epsilon())
250  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
251 
252  // Calculate displacement in x from relative displacement in y-z plane
253  const CartesianVector firstPositionYZ(0.f, firstVertex.GetPosition().GetY(), firstVertex.GetPosition().GetZ());
254  const CartesianVector firstDirectionYZ(0.f, firstVertex.GetDirection().GetY(), firstVertex.GetDirection().GetZ());
255 
256  const CartesianVector secondPositionYZ(0.f, secondVertex.GetPosition().GetY(), secondVertex.GetPosition().GetZ());
257  const CartesianVector secondDirectionYZ(0.f, secondVertex.GetDirection().GetY(), secondVertex.GetDirection().GetZ());
258 
259  const float firstDirectionYZmag(firstDirectionYZ.GetMagnitude());
260  const float secondDirectionYZmag(secondDirectionYZ.GetMagnitude());
261 
262  if (firstDirectionYZmag < std::numeric_limits<float>::epsilon() || secondDirectionYZmag < std::numeric_limits<float>::epsilon())
263  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
264 
265  const float R1(firstDirectionYZ.GetUnitVector().GetDotProduct(firstPositionYZ - secondPositionYZ) / firstDirectionYZmag);
266  const float X1(-1.f *
267  firstDirectionX.GetUnitVector().GetDotProduct(secondVertex.GetPosition() - (firstVertex.GetPosition() - firstVertex.GetDirection() * R1)));
268 
269  const float R2(secondDirectionYZ.GetUnitVector().GetDotProduct(secondPositionYZ - firstPositionYZ) / secondDirectionYZmag);
270  const float X2(-1.f *
271  secondDirectionX.GetUnitVector().GetDotProduct(firstVertex.GetPosition() - (secondVertex.GetPosition() - secondVertex.GetDirection() * R2)));
272 
273  // ATTN: By convention, X0 is half the displacement in x (because both Pfos will be corrected)
274  return (X1 + X2) * 0.25f;
275 }
276 
277 //------------------------------------------------------------------------------------------------------------------------------------------
278 
279 bool LArStitchingHelper::SortTPCs(const pandora::LArTPC *const pLhs, const pandora::LArTPC *const pRhs)
280 {
281  if (std::fabs(pLhs->GetCenterX() - pRhs->GetCenterX()) > std::numeric_limits<float>::epsilon())
282  return (pLhs->GetCenterX() < pRhs->GetCenterX());
283 
284  if (std::fabs(pLhs->GetCenterY() - pRhs->GetCenterY()) > std::numeric_limits<float>::epsilon())
285  return (pLhs->GetCenterY() < pRhs->GetCenterY());
286 
287  return (pLhs->GetCenterZ() < pRhs->GetCenterZ());
288 }
289 
290 //------------------------------------------------------------------------------------------------------------------------------------------
291 
292 bool LArStitchingHelper::HasPfoBeenStitched(const ParticleFlowObject *const pPfo)
293 {
294  const PropertiesMap &properties(pPfo->GetPropertiesMap());
295  const PropertiesMap::const_iterator iter(properties.find("X0"));
296 
297  if (iter != properties.end())
298  return true;
299 
300  return false;
301 }
302 
303 //------------------------------------------------------------------------------------------------------------------------------------------
304 
305 float LArStitchingHelper::GetPfoX0(const ParticleFlowObject *const pPfo)
306 {
307  // ATTN: If no stitch is present return a shift of 0.f
308  if (!LArStitchingHelper::HasPfoBeenStitched(pPfo))
309  return 0.f;
310 
311  // ATTN: HasPfoBeenStitched ensures X0 is present in the properties map
312  return pPfo->GetPropertiesMap().at("X0");
313 }
314 
315 } // namespace lar_content
intermediate_table::const_iterator const_iterator
LArPointingCluster class.
Double_t X2
Definition: plot.C:263
TFile f
Definition: plotHisto.C:6
Double_t X1
Definition: plot.C:263
const Vertex & GetOuterVertex() const
Get the outer vertex.
const Vertex & GetInnerVertex() const
Get the inner vertex.
Header file for the MultiPandoraApi class.
const pandora::CartesianVector & GetDirection() const
Get the vertex direction.
Header file for the helper class for multiple drift volumes.
static const PandoraInstanceMap & GetPandoraInstanceMap()
Get the pandora instance map.
const pandora::CartesianVector & GetPosition() const
Get the vertex position.