LArSoft  v07_13_02
Liquid Argon Software toolkit - http://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,
179  const LArPointingCluster &pointingCluster1, const LArPointingCluster &pointingCluster2,
180  LArPointingCluster::Vertex &closestVertex1, LArPointingCluster::Vertex &closestVertex2)
181 {
182  if (&larTPC1 == &larTPC2)
183  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
184 
185  // Find the closest vertices based on relative X positions in tpc
186  const float dxVolume(larTPC2.GetCenterX() - larTPC1.GetCenterX());
187  const float dx1(pointingCluster1.GetOuterVertex().GetPosition().GetX() - pointingCluster1.GetInnerVertex().GetPosition().GetX());
188  const float dx2(pointingCluster2.GetOuterVertex().GetPosition().GetX() - pointingCluster2.GetInnerVertex().GetPosition().GetX());
189 
190  if (std::fabs(dx1) < std::numeric_limits<float>::epsilon() || std::fabs(dx2) < std::numeric_limits<float>::epsilon())
191  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
192 
193  const bool useInner1((dxVolume > 0.f) == (dx1 < 0.f)); // xVol1 - OUTER - INNER - | - xVol2 [xVol2-xVol1>0; xOuter1-xInner1<0]
194  const bool useInner2((dxVolume > 0.f) == (dx2 > 0.f)); // xVol1 - | - INNER - OUTER - xVol2 [xVol2-xVol1>0; xOuter2-xInner2>0]
195 
196  // Confirm that these really are the closest pair of vertices by checking the other possible pairs
197  const LArPointingCluster::Vertex &nearVertex1(useInner1 ? pointingCluster1.GetInnerVertex() : pointingCluster1.GetOuterVertex());
198  const LArPointingCluster::Vertex &nearVertex2(useInner2 ? pointingCluster2.GetInnerVertex() : pointingCluster2.GetOuterVertex());
199 
200  const LArPointingCluster::Vertex &farVertex1(useInner1 ? pointingCluster1.GetOuterVertex() : pointingCluster1.GetInnerVertex());
201  const LArPointingCluster::Vertex &farVertex2(useInner2 ? pointingCluster2.GetOuterVertex() : pointingCluster2.GetInnerVertex());
202 
203  const float dxNearNear(0.f);
204  const float dyNearNear(nearVertex1.GetPosition().GetY() - nearVertex2.GetPosition().GetY());
205  const float dzNearNear(nearVertex1.GetPosition().GetZ() - nearVertex2.GetPosition().GetZ());
206  const float drNearNearSquared(dxNearNear * dxNearNear + dyNearNear * dyNearNear + dzNearNear * dzNearNear);
207 
208  const float dxNearFar(std::fabs(dx2));
209  const float dyNearFar(nearVertex1.GetPosition().GetY() - farVertex2.GetPosition().GetY());
210  const float dzNearFar(nearVertex1.GetPosition().GetZ() - farVertex2.GetPosition().GetZ());
211  const float drNearFarSquared(dxNearFar * dxNearFar + dyNearFar * dyNearFar + dzNearFar * dzNearFar);
212 
213  const float dxFarNear(std::fabs(dx1));
214  const float dyFarNear(farVertex1.GetPosition().GetY() - nearVertex2.GetPosition().GetY());
215  const float dzFarNear(farVertex1.GetPosition().GetZ() - nearVertex2.GetPosition().GetZ());
216  const float drFarNearSquared(dxFarNear * dxFarNear + dyFarNear * dyFarNear + dzFarNear * dzFarNear);
217 
218  const float dxFarFar(std::fabs(dx1) + std::fabs(dx2));
219  const float dyFarFar(farVertex1.GetPosition().GetY() - farVertex2.GetPosition().GetY());
220  const float dzFarFar(farVertex1.GetPosition().GetZ() - farVertex2.GetPosition().GetZ());
221  const float drFarFarSquared(dxFarFar * dxFarFar + dyFarFar * dyFarFar + dzFarFar * dzFarFar);
222 
223  if (drNearNearSquared > std::min(drFarFarSquared, std::min(drNearFarSquared, drFarNearSquared)))
224  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
225 
226  closestVertex1 = nearVertex1;
227  closestVertex2 = nearVertex2;
228 }
229 
230 //------------------------------------------------------------------------------------------------------------------------------------------
231 
232 float LArStitchingHelper::CalculateX0(const LArTPC &firstTPC, const LArTPC &secondTPC,
233  const LArPointingCluster::Vertex &firstVertex, const LArPointingCluster::Vertex &secondVertex)
234 {
235  if (&firstTPC == &secondTPC)
236  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
237 
238  // ATTN: Assume that Pfos are crossing either an anode-anode boundary or a cathode-cathode boundary
239  if (firstTPC.IsDriftInPositiveX() == secondTPC.IsDriftInPositiveX())
240  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
241 
242  // Assume that Pfos have opposite direction components in x, and have some direction component in the y-z plane
243  const CartesianVector firstDirectionX(firstVertex.GetDirection().GetX(), 0.f, 0.f);
244  const CartesianVector secondDirectionX(secondVertex.GetDirection().GetX(), 0.f, 0.f);
245 
246  if (std::fabs(firstDirectionX.GetX()) > 1.f - std::numeric_limits<float>::epsilon() ||
247  std::fabs(secondDirectionX.GetX()) > 1.f - std::numeric_limits<float>::epsilon())
248  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
249 
250  if (-firstDirectionX.GetDotProduct(secondDirectionX) < std::numeric_limits<float>::epsilon())
251  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
252 
253  // Calculate displacement in x from relative displacement in y-z plane
254  const CartesianVector firstPositionYZ(0.f, firstVertex.GetPosition().GetY(), firstVertex.GetPosition().GetZ());
255  const CartesianVector firstDirectionYZ(0.f, firstVertex.GetDirection().GetY(), firstVertex.GetDirection().GetZ());
256 
257  const CartesianVector secondPositionYZ(0.f, secondVertex.GetPosition().GetY(), secondVertex.GetPosition().GetZ());
258  const CartesianVector secondDirectionYZ(0.f, secondVertex.GetDirection().GetY(), secondVertex.GetDirection().GetZ());
259 
260  const float firstDirectionYZmag(firstDirectionYZ.GetMagnitude());
261  const float secondDirectionYZmag(secondDirectionYZ.GetMagnitude());
262 
263  if (firstDirectionYZmag < std::numeric_limits<float>::epsilon() || secondDirectionYZmag < std::numeric_limits<float>::epsilon())
264  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
265 
266  const float R1(firstDirectionYZ.GetUnitVector().GetDotProduct(firstPositionYZ - secondPositionYZ) / firstDirectionYZmag);
267  const float X1(-1.f * 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 * secondDirectionX.GetUnitVector().GetDotProduct(firstVertex.GetPosition() - (secondVertex.GetPosition() - secondVertex.GetDirection() * R2)));
271 
272  // ATTN: By convention, X0 is half the displacement in x (because both Pfos will be corrected)
273  return (X1 + X2) * 0.25f;
274 }
275 
276 //------------------------------------------------------------------------------------------------------------------------------------------
277 
278 CartesianVector LArStitchingHelper::GetCorrectedPosition(const LArTPC &larTPC, const float x0, const CartesianVector &inputPosition)
279 {
280  return (inputPosition + CartesianVector(larTPC.IsDriftInPositiveX() ? -x0 : x0, 0.f, 0.f));
281 }
282 
283 //------------------------------------------------------------------------------------------------------------------------------------------
284 
285 bool LArStitchingHelper::SortTPCs(const pandora::LArTPC *const pLhs, const pandora::LArTPC *const pRhs)
286 {
287  if (std::fabs(pLhs->GetCenterX() - pRhs->GetCenterX()) > std::numeric_limits<float>::epsilon())
288  return (pLhs->GetCenterX() < pRhs->GetCenterX());
289 
290  if (std::fabs(pLhs->GetCenterY() - pRhs->GetCenterY()) > std::numeric_limits<float>::epsilon())
291  return (pLhs->GetCenterY() < pRhs->GetCenterY());
292 
293  return (pLhs->GetCenterZ() < pRhs->GetCenterZ());
294 }
295 
296 } // namespace lar_content
LArPointingCluster class.
Double_t X2
Definition: plot.C:266
TFile f
Definition: plotHisto.C:6
Int_t max
Definition: plot.C:27
Double_t X1
Definition: plot.C:266
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.
Int_t min
Definition: plot.C:26
static const PandoraInstanceMap & GetPandoraInstanceMap()
Get the pandora instance map.
const pandora::CartesianVector & GetPosition() const
Get the vertex position.