LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArGeometryHelper.cc
Go to the documentation of this file.
1 
9 #include "Managers/GeometryManager.h"
10 #include "Managers/PluginManager.h"
11 
12 #include "Geometry/DetectorGap.h"
13 #include "Geometry/LArTPC.h"
14 
15 #include "Objects/CartesianVector.h"
16 
17 #include "Pandora/Pandora.h"
18 
21 
23 
24 #include "Plugins/LArTransformationPlugin.h"
25 
26 using namespace pandora;
27 
28 namespace lar_content
29 {
30 
31 float LArGeometryHelper::MergeTwoPositions(const Pandora &pandora, const HitType view1, const HitType view2, const float position1, const float position2)
32 {
33  if (view1 == view2)
34  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
35 
36  if ((view1 == TPC_VIEW_U) && (view2 == TPC_VIEW_V))
37  {
38  return pandora.GetPlugins()->GetLArTransformationPlugin()->UVtoW(position1, position2);
39  }
40 
41  if ((view1 == TPC_VIEW_V) && (view2 == TPC_VIEW_U))
42  {
43  return pandora.GetPlugins()->GetLArTransformationPlugin()->UVtoW(position2, position1);
44  }
45 
46  if ((view1 == TPC_VIEW_W) && (view2 == TPC_VIEW_U))
47  {
48  return pandora.GetPlugins()->GetLArTransformationPlugin()->WUtoV(position1, position2);
49  }
50 
51  if ((view1 == TPC_VIEW_U) && (view2 == TPC_VIEW_W))
52  {
53  return pandora.GetPlugins()->GetLArTransformationPlugin()->WUtoV(position2, position1);
54  }
55 
56  if ((view1 == TPC_VIEW_V) && (view2 == TPC_VIEW_W))
57  {
58  return pandora.GetPlugins()->GetLArTransformationPlugin()->VWtoU(position1, position2);
59  }
60 
61  if ((view1 == TPC_VIEW_W) && (view2 == TPC_VIEW_V))
62  {
63  return pandora.GetPlugins()->GetLArTransformationPlugin()->VWtoU(position2, position1);
64  }
65 
66  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
67 }
68 
69 //------------------------------------------------------------------------------------------------------------------------------------------
70 
71 CartesianVector LArGeometryHelper::MergeTwoDirections(
72  const Pandora &pandora, const HitType view1, const HitType view2, const CartesianVector &direction1, const CartesianVector &direction2)
73 {
74  // x components must be consistent
75  if (direction1.GetX() * direction2.GetX() < 0.f)
76  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
77 
78  // scale x components to a common value
79  const float s1((std::fabs(direction1.GetX()) > std::numeric_limits<float>::epsilon()) ? 100.f * std::fabs(direction2.GetX()) : 1.f);
80  const float s2((std::fabs(direction2.GetX()) > std::numeric_limits<float>::epsilon()) ? 100.f * std::fabs(direction1.GetX()) : 1.f);
81 
82  float pX(s1 * direction1.GetX()), pU(0.f), pV(0.f), pW(0.f);
83  HitType newView(HIT_CUSTOM);
84 
85  if ((view1 == TPC_VIEW_U) && (view2 == TPC_VIEW_V))
86  {
87  pU = s1 * direction1.GetZ();
88  pV = s2 * direction2.GetZ();
89  newView = TPC_VIEW_W;
90  }
91 
92  if ((view1 == TPC_VIEW_V) && (view2 == TPC_VIEW_U))
93  {
94  pV = s1 * direction1.GetZ();
95  pU = s2 * direction2.GetZ();
96  newView = TPC_VIEW_W;
97  }
98 
99  if ((view1 == TPC_VIEW_W) && (view2 == TPC_VIEW_U))
100  {
101  pW = s1 * direction1.GetZ();
102  pU = s2 * direction2.GetZ();
103  newView = TPC_VIEW_V;
104  }
105 
106  if ((view1 == TPC_VIEW_U) && (view2 == TPC_VIEW_W))
107  {
108  pU = s1 * direction1.GetZ();
109  pW = s2 * direction2.GetZ();
110  newView = TPC_VIEW_V;
111  }
112 
113  if ((view1 == TPC_VIEW_V) && (view2 == TPC_VIEW_W))
114  {
115  pV = s1 * direction1.GetZ();
116  pW = s2 * direction2.GetZ();
117  newView = TPC_VIEW_U;
118  }
119 
120  if ((view1 == TPC_VIEW_W) && (view2 == TPC_VIEW_V))
121  {
122  pW = s1 * direction1.GetZ();
123  pV = s2 * direction2.GetZ();
124  newView = TPC_VIEW_U;
125  }
126 
127  if (newView == TPC_VIEW_W)
128  return CartesianVector(pX, 0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->UVtoW(pU, pV)).GetUnitVector();
129 
130  if (newView == TPC_VIEW_U)
131  return CartesianVector(pX, 0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->VWtoU(pV, pW)).GetUnitVector();
132 
133  if (newView == TPC_VIEW_V)
134  return CartesianVector(pX, 0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->WUtoV(pW, pU)).GetUnitVector();
135 
136  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
137 }
138 
139 //------------------------------------------------------------------------------------------------------------------------------------------
140 
141 void LArGeometryHelper::MergeTwoPositions(const Pandora &pandora, const HitType view1, const HitType view2,
142  const CartesianVector &position1, const CartesianVector &position2, CartesianVector &position3, float &chiSquared)
143 {
144  if (view1 == view2)
145  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
146 
147  const float X3((position1.GetX() + position2.GetX()) / 2.f);
148  const float Z1(position1.GetZ());
149  const float Z2(position2.GetZ());
150  const float Z3(LArGeometryHelper::MergeTwoPositions(pandora, view1, view2, Z1, Z2));
151 
152  position3.SetValues(X3, 0.f, Z3);
153  const float sigmaUVW(LArGeometryHelper::GetSigmaUVW(pandora));
154  chiSquared = ((X3 - position1.GetX()) * (X3 - position1.GetX()) + (X3 - position2.GetX()) * (X3 - position2.GetX())) / (sigmaUVW * sigmaUVW);
155 }
156 
157 //------------------------------------------------------------------------------------------------------------------------------------------
158 
159 void LArGeometryHelper::MergeTwoPositions(const Pandora &pandora, const HitType view1, const HitType view2, const CartesianVector &position1,
160  const CartesianVector &position2, CartesianVector &outputU, CartesianVector &outputV, CartesianVector &outputW, float &chiSquared)
161 {
162  if (view1 == view2)
163  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
164 
165  CartesianVector position3(0.f, 0.f, 0.f);
166  LArGeometryHelper::MergeTwoPositions(pandora, view1, view2, position1, position2, position3, chiSquared);
167 
168  float aveU(0.f), aveV(0.f), aveW(0.f);
169  const float Z1(position1.GetZ()), Z2(position2.GetZ()), Z3(position3.GetZ()), aveX(position3.GetX());
170 
171  if ((view1 == TPC_VIEW_U) && (view2 == TPC_VIEW_V))
172  {
173  aveU = Z1;
174  aveV = Z2;
175  aveW = Z3;
176  }
177 
178  if ((view1 == TPC_VIEW_V) && (view2 == TPC_VIEW_W))
179  {
180  aveV = Z1;
181  aveW = Z2;
182  aveU = Z3;
183  }
184 
185  if ((view1 == TPC_VIEW_W) && (view2 == TPC_VIEW_U))
186  {
187  aveW = Z1;
188  aveU = Z2;
189  aveV = Z3;
190  }
191 
192  if ((view1 == TPC_VIEW_V) && (view2 == TPC_VIEW_U))
193  {
194  aveV = Z1;
195  aveU = Z2;
196  aveW = Z3;
197  }
198 
199  if ((view1 == TPC_VIEW_W) && (view2 == TPC_VIEW_V))
200  {
201  aveW = Z1;
202  aveV = Z2;
203  aveU = Z3;
204  }
205 
206  if ((view1 == TPC_VIEW_U) && (view2 == TPC_VIEW_W))
207  {
208  aveU = Z1;
209  aveW = Z2;
210  aveV = Z3;
211  }
212 
213  outputU.SetValues(aveX, 0.f, aveU);
214  outputV.SetValues(aveX, 0.f, aveV);
215  outputW.SetValues(aveX, 0.f, aveW);
216 }
217 
218 //------------------------------------------------------------------------------------------------------------------------------------------
219 
220 void LArGeometryHelper::MergeThreePositions(const Pandora &pandora, const HitType view1, const HitType view2, const HitType view3,
221  const CartesianVector &position1, const CartesianVector &position2, const CartesianVector &position3, CartesianVector &outputU,
222  CartesianVector &outputV, CartesianVector &outputW, float &chiSquared)
223 {
224  if ((view1 == view2) || (view2 == view3) || (view3 == view1))
225  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
226 
227  if ((view1 == TPC_VIEW_U) && (view2 == TPC_VIEW_V))
228  {
229  return LArGeometryHelper::MergeThreePositions(pandora, position1, position2, position3, outputU, outputV, outputW, chiSquared);
230  }
231 
232  if ((view1 == TPC_VIEW_V) && (view2 == TPC_VIEW_W))
233  {
234  return LArGeometryHelper::MergeThreePositions(pandora, position3, position1, position2, outputU, outputV, outputW, chiSquared);
235  }
236 
237  if ((view1 == TPC_VIEW_W) && (view2 == TPC_VIEW_U))
238  {
239  return LArGeometryHelper::MergeThreePositions(pandora, position2, position3, position1, outputU, outputV, outputW, chiSquared);
240  }
241 
242  if ((view1 == TPC_VIEW_V) && (view2 == TPC_VIEW_U))
243  {
244  return LArGeometryHelper::MergeThreePositions(pandora, position2, position1, position3, outputU, outputV, outputW, chiSquared);
245  }
246 
247  if ((view1 == TPC_VIEW_W) && (view2 == TPC_VIEW_V))
248  {
249  return LArGeometryHelper::MergeThreePositions(pandora, position3, position2, position1, outputU, outputV, outputW, chiSquared);
250  }
251 
252  if ((view1 == TPC_VIEW_U) && (view2 == TPC_VIEW_W))
253  {
254  return LArGeometryHelper::MergeThreePositions(pandora, position1, position3, position2, outputU, outputV, outputW, chiSquared);
255  }
256 
257  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
258 }
259 
260 //------------------------------------------------------------------------------------------------------------------------------------------
261 
262 void LArGeometryHelper::MergeThreePositions(const Pandora &pandora, const CartesianVector &positionU, const CartesianVector &positionV,
263  const CartesianVector &positionW, CartesianVector &outputU, CartesianVector &outputV, CartesianVector &outputW, float &chiSquared)
264 {
265  const float YfromUV(pandora.GetPlugins()->GetLArTransformationPlugin()->UVtoY(positionU.GetZ(), positionV.GetZ()));
266  const float YfromUW(pandora.GetPlugins()->GetLArTransformationPlugin()->UWtoY(positionU.GetZ(), positionW.GetZ()));
267  const float YfromVW(pandora.GetPlugins()->GetLArTransformationPlugin()->VWtoY(positionV.GetZ(), positionW.GetZ()));
268 
269  const float ZfromUV(pandora.GetPlugins()->GetLArTransformationPlugin()->UVtoZ(positionU.GetZ(), positionV.GetZ()));
270  const float ZfromUW(pandora.GetPlugins()->GetLArTransformationPlugin()->UWtoZ(positionU.GetZ(), positionW.GetZ()));
271  const float ZfromVW(pandora.GetPlugins()->GetLArTransformationPlugin()->VWtoZ(positionV.GetZ(), positionW.GetZ()));
272 
273  // ATTN For detectors where w and z are equivalent, remain consistent with original treatment. TODO Use new treatment always.
274  const bool useOldWZEquivalentTreatment(std::fabs(ZfromUW - ZfromVW) < std::numeric_limits<float>::epsilon());
275  const float aveX((positionU.GetX() + positionV.GetX() + positionW.GetX()) / 3.f);
276  const float aveY(useOldWZEquivalentTreatment ? YfromUV : (YfromUV + YfromUW + YfromVW) / 3.f);
277  const float aveZ(useOldWZEquivalentTreatment ? (positionW.GetZ() + 2.f * ZfromUV) / 3.f : (ZfromUV + ZfromUW + ZfromVW) / 3.f);
278 
279  const float aveU(pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoU(aveY, aveZ));
280  const float aveV(pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoV(aveY, aveZ));
281  const float aveW(pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoW(aveY, aveZ));
282 
283  outputU.SetValues(aveX, 0.f, aveU);
284  outputV.SetValues(aveX, 0.f, aveV);
285  outputW.SetValues(aveX, 0.f, aveW);
286 
287  const float sigmaUVW(LArGeometryHelper::GetSigmaUVW(pandora));
288  chiSquared = ((outputU.GetX() - positionU.GetX()) * (outputU.GetX() - positionU.GetX()) +
289  (outputV.GetX() - positionV.GetX()) * (outputV.GetX() - positionV.GetX()) +
290  (outputW.GetX() - positionW.GetX()) * (outputW.GetX() - positionW.GetX()) +
291  (outputU.GetZ() - positionU.GetZ()) * (outputU.GetZ() - positionU.GetZ()) +
292  (outputV.GetZ() - positionV.GetZ()) * (outputV.GetZ() - positionV.GetZ()) +
293  (outputW.GetZ() - positionW.GetZ()) * (outputW.GetZ() - positionW.GetZ())) /
294  (sigmaUVW * sigmaUVW);
295 }
296 
297 //------------------------------------------------------------------------------------------------------------------------------------------
298 
299 void LArGeometryHelper::MergeTwoPositions3D(const Pandora &pandora, const HitType view1, const HitType view2,
300  const CartesianVector &position1, const CartesianVector &position2, CartesianVector &position3D, float &chiSquared)
301 {
302  CartesianVector positionU(0.f, 0.f, 0.f), positionV(0.f, 0.f, 0.f), positionW(0.f, 0.f, 0.f);
303  LArGeometryHelper::MergeTwoPositions(pandora, view1, view2, position1, position2, positionU, positionV, positionW, chiSquared);
304 
305  position3D.SetValues(positionW.GetX(), pandora.GetPlugins()->GetLArTransformationPlugin()->UVtoY(positionU.GetZ(), positionV.GetZ()),
306  pandora.GetPlugins()->GetLArTransformationPlugin()->UVtoZ(positionU.GetZ(), positionV.GetZ()));
307 }
308 
309 //------------------------------------------------------------------------------------------------------------------------------------------
310 
311 void LArGeometryHelper::MergeThreePositions3D(const Pandora &pandora, const HitType view1, const HitType view2, const HitType view3,
312  const CartesianVector &position1, const CartesianVector &position2, const CartesianVector &position3, CartesianVector &position3D, float &chiSquared)
313 {
314  CartesianVector positionU(0.f, 0.f, 0.f), positionV(0.f, 0.f, 0.f), positionW(0.f, 0.f, 0.f);
315  LArGeometryHelper::MergeThreePositions(pandora, view1, view2, view3, position1, position2, position3, positionU, positionV, positionW, chiSquared);
316 
317  position3D.SetValues(positionW.GetX(), pandora.GetPlugins()->GetLArTransformationPlugin()->UVtoY(positionU.GetZ(), positionV.GetZ()),
318  pandora.GetPlugins()->GetLArTransformationPlugin()->UVtoZ(positionU.GetZ(), positionV.GetZ()));
319 }
320 
321 //------------------------------------------------------------------------------------------------------------------------------------------
322 
323 CartesianVector LArGeometryHelper::ProjectPosition(const Pandora &pandora, const CartesianVector &position3D, const HitType view)
324 {
325  if (view == TPC_VIEW_U)
326  {
327  return CartesianVector(
328  position3D.GetX(), 0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoU(position3D.GetY(), position3D.GetZ()));
329  }
330 
331  else if (view == TPC_VIEW_V)
332  {
333  return CartesianVector(
334  position3D.GetX(), 0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoV(position3D.GetY(), position3D.GetZ()));
335  }
336 
337  else if (view == TPC_VIEW_W)
338  {
339  return CartesianVector(
340  position3D.GetX(), 0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoW(position3D.GetY(), position3D.GetZ()));
341  }
342 
343  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
344 }
345 
346 //------------------------------------------------------------------------------------------------------------------------------------------
347 
348 CartesianVector LArGeometryHelper::ProjectDirection(const Pandora &pandora, const CartesianVector &direction3D, const HitType view)
349 {
350  if (view == TPC_VIEW_U)
351  {
352  return CartesianVector(
353  direction3D.GetX(), 0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoU(direction3D.GetY(), direction3D.GetZ()))
354  .GetUnitVector();
355  }
356 
357  else if (view == TPC_VIEW_V)
358  {
359  return CartesianVector(
360  direction3D.GetX(), 0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoV(direction3D.GetY(), direction3D.GetZ()))
361  .GetUnitVector();
362  }
363 
364  else if (view == TPC_VIEW_W)
365  {
366  return CartesianVector(
367  direction3D.GetX(), 0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoW(direction3D.GetY(), direction3D.GetZ()))
368  .GetUnitVector();
369  }
370 
371  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
372 }
373 
374 //------------------------------------------------------------------------------------------------------------------------------------------
375 
376 float LArGeometryHelper::GetWirePitch(const Pandora &pandora, const HitType view, const float maxWirePitchDiscrepancy)
377 {
378  if (view != TPC_VIEW_U && view != TPC_VIEW_V && view != TPC_VIEW_W && view != TPC_3D)
379  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
380 
381  if (view == TPC_3D)
382  {
383  const float pitchU{LArGeometryHelper::GetWirePitch(pandora, TPC_VIEW_U, maxWirePitchDiscrepancy)};
384  const float pitchV{LArGeometryHelper::GetWirePitch(pandora, TPC_VIEW_V, maxWirePitchDiscrepancy)};
385  const float pitchW{LArGeometryHelper::GetWirePitch(pandora, TPC_VIEW_W, maxWirePitchDiscrepancy)};
386  return std::max({pitchU, pitchV, pitchW});
387  }
388 
389  const LArTPCMap &larTPCMap(pandora.GetGeometry()->GetLArTPCMap());
390 
391  if (larTPCMap.empty())
392  {
393  std::cout << "LArGeometryHelper::GetWirePitch - LArTPC description not registered with Pandora as required " << std::endl;
394  throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
395  }
396 
397  const LArTPC *const pFirstLArTPC(larTPCMap.begin()->second);
398  const float wirePitch(view == TPC_VIEW_U ? pFirstLArTPC->GetWirePitchU()
399  : (view == TPC_VIEW_V ? pFirstLArTPC->GetWirePitchV() : pFirstLArTPC->GetWirePitchW()));
400 
401  for (const LArTPCMap::value_type &mapEntry : larTPCMap)
402  {
403  const LArTPC *const pLArTPC(mapEntry.second);
404  const float alternateWirePitch(
405  view == TPC_VIEW_U ? pLArTPC->GetWirePitchU() : (view == TPC_VIEW_V ? pLArTPC->GetWirePitchV() : pLArTPC->GetWirePitchW()));
406 
407  if (std::fabs(wirePitch - alternateWirePitch) > maxWirePitchDiscrepancy)
408  {
409  std::cout << "LArGeometryHelper::GetWirePitch - LArTPC configuration not supported" << std::endl;
410  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
411  }
412  }
413 
414  return wirePitch;
415 }
416 
417 //------------------------------------------------------------------------------------------------------------------------------------------
418 
419 float LArGeometryHelper::GetWirePitchRatio(const Pandora &pandora, const HitType view)
420 {
421  const float pitchU{LArGeometryHelper::GetWirePitch(pandora, TPC_VIEW_U)};
422  const float pitchV{LArGeometryHelper::GetWirePitch(pandora, TPC_VIEW_V)};
423  const float pitchW{LArGeometryHelper::GetWirePitch(pandora, TPC_VIEW_W)};
424  const float pitchMin{std::min({pitchU, pitchV, pitchW})};
425  const float pitchView{LArGeometryHelper::GetWirePitch(pandora, view)};
426 
427  return pitchView / pitchMin;
428 }
429 
430 //------------------------------------------------------------------------------------------------------------------------------------------
431 
432 CartesianVector LArGeometryHelper::GetWireAxis(const Pandora &pandora, const HitType view)
433 {
434  if (view == TPC_VIEW_U)
435  {
436  return CartesianVector(0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoU(1.f, 0.f),
437  pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoU(0.f, 1.f));
438  }
439 
440  else if (view == TPC_VIEW_V)
441  {
442  return CartesianVector(0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoV(1.f, 0.f),
443  pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoV(0.f, 1.f));
444  }
445 
446  else if (view == TPC_VIEW_W)
447  {
448  return CartesianVector(0.f, pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoW(1.f, 0.f),
449  pandora.GetPlugins()->GetLArTransformationPlugin()->YZtoW(0.f, 1.f));
450  }
451 
452  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
453 }
454 
455 //------------------------------------------------------------------------------------------------------------------------------------------
456 
457 void LArGeometryHelper::GetCommonDaughterVolumes(const Cluster *const pCluster1, const Cluster *const pCluster2, UIntSet &intersect)
458 {
459  UIntSet daughterVolumeIds1, daughterVolumeIds2;
460 
461  LArClusterHelper::GetDaughterVolumeIDs(pCluster1, daughterVolumeIds1);
462  LArClusterHelper::GetDaughterVolumeIDs(pCluster2, daughterVolumeIds2);
463 
464  std::set_intersection(daughterVolumeIds1.begin(), daughterVolumeIds1.end(), daughterVolumeIds2.begin(), daughterVolumeIds2.end(),
465  std::inserter(intersect, intersect.begin()));
466 }
467 
468 //------------------------------------------------------------------------------------------------------------------------------------------
469 
470 bool LArGeometryHelper::IsInGap(const Pandora &pandora, const CartesianVector &testPoint2D, const HitType hitType, const float gapTolerance)
471 {
472  // ATTN: input test point MUST be a 2D position vector
473  for (const DetectorGap *const pDetectorGap : pandora.GetGeometry()->GetDetectorGapList())
474  {
475  if (pDetectorGap->IsInGap(testPoint2D, hitType, gapTolerance))
476  return true;
477  }
478 
479  return false;
480 }
481 
482 //------------------------------------------------------------------------------------------------------------------------------------------
483 
484 bool LArGeometryHelper::IsInGap3D(const Pandora &pandora, const CartesianVector &testPoint3D, const HitType hitType, const float gapTolerance)
485 {
486  const CartesianVector testPoint2D(LArGeometryHelper::ProjectPosition(pandora, testPoint3D, hitType));
487  return LArGeometryHelper::IsInGap(pandora, testPoint2D, hitType, gapTolerance);
488 }
489 
490 //------------------------------------------------------------------------------------------------------------------------------------------
491 
492 bool LArGeometryHelper::IsXSamplingPointInGap(const Pandora &pandora, const float xSample, const TwoDSlidingFitResult &slidingFitResult, const float gapTolerance)
493 {
494  const HitType hitType(LArClusterHelper::GetClusterHitType(slidingFitResult.GetCluster()));
495  const CartesianVector minLayerPosition(slidingFitResult.GetGlobalMinLayerPosition());
496  const CartesianVector maxLayerPosition(slidingFitResult.GetGlobalMaxLayerPosition());
497 
498  const bool minLayerIsAtLowX(minLayerPosition.GetX() < maxLayerPosition.GetX());
499  const CartesianVector &lowXCoordinate(minLayerIsAtLowX ? minLayerPosition : maxLayerPosition);
500  const CartesianVector &highXCoordinate(minLayerIsAtLowX ? maxLayerPosition : minLayerPosition);
501 
502  if ((xSample > lowXCoordinate.GetX()) && (xSample < highXCoordinate.GetX()))
503  {
504  CartesianVector slidingFitPosition(0.f, 0.f, 0.f);
505 
506  if (STATUS_CODE_SUCCESS == slidingFitResult.GetGlobalFitPositionAtX(xSample, slidingFitPosition))
507  return (LArGeometryHelper::IsInGap(pandora, slidingFitPosition, hitType, gapTolerance));
508  }
509 
510  const CartesianVector lowXDirection(
511  minLayerIsAtLowX ? slidingFitResult.GetGlobalMinLayerDirection() : slidingFitResult.GetGlobalMaxLayerDirection());
512  const CartesianVector highXDirection(
513  minLayerIsAtLowX ? slidingFitResult.GetGlobalMaxLayerDirection() : slidingFitResult.GetGlobalMinLayerDirection());
514 
515  const bool sampleIsNearerToLowX(std::fabs(xSample - lowXCoordinate.GetX()) < std::fabs(xSample - highXCoordinate.GetX()));
516  const CartesianVector &startPosition(sampleIsNearerToLowX ? lowXCoordinate : highXCoordinate);
517  const CartesianVector &startDirection(sampleIsNearerToLowX ? lowXDirection : highXDirection);
518 
519  if (std::fabs(startDirection.GetX()) < std::numeric_limits<float>::epsilon())
520  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
521 
522  const float pathLength((xSample - startPosition.GetX()) / startDirection.GetX());
523  const CartesianVector samplingPoint(startPosition + startDirection * pathLength);
524 
525  return (LArGeometryHelper::IsInGap(pandora, samplingPoint, hitType, gapTolerance));
526 }
527 
528 //------------------------------------------------------------------------------------------------------------------------------------------
529 
530 float LArGeometryHelper::CalculateGapDeltaZ(const Pandora &pandora, const float minZ, const float maxZ, const HitType hitType)
531 {
532  if (maxZ - minZ < std::numeric_limits<float>::epsilon())
533  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
534 
535  float gapDeltaZ(0.f);
536 
537  for (const DetectorGap *const pDetectorGap : pandora.GetGeometry()->GetDetectorGapList())
538  {
539  const LineGap *const pLineGap = dynamic_cast<const LineGap *>(pDetectorGap);
540 
541  if (!pLineGap)
542  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
543 
544  const LineGapType lineGapType(pLineGap->GetLineGapType());
545 
546  if (!(((TPC_VIEW_U == hitType) && (TPC_WIRE_GAP_VIEW_U == lineGapType)) || ((TPC_VIEW_V == hitType) && (TPC_WIRE_GAP_VIEW_V == lineGapType)) ||
547  ((TPC_VIEW_W == hitType) && (TPC_WIRE_GAP_VIEW_W == lineGapType))))
548  {
549  continue;
550  }
551 
552  if ((pLineGap->GetLineStartZ() > maxZ) || (pLineGap->GetLineEndZ() < minZ))
553  continue;
554 
555  const float gapMinZ(std::max(minZ, pLineGap->GetLineStartZ()));
556  const float gapMaxZ(std::min(maxZ, pLineGap->GetLineEndZ()));
557 
558  if ((gapMaxZ - gapMinZ) > std::numeric_limits<float>::epsilon())
559  gapDeltaZ += (gapMaxZ - gapMinZ);
560  }
561 
562  return gapDeltaZ;
563 }
564 
565 //------------------------------------------------------------------------------------------------------------------------------------------
566 
567 float LArGeometryHelper::GetSigmaUVW(const Pandora &pandora, const float maxSigmaDiscrepancy)
568 {
569  const LArTPCMap &larTPCMap(pandora.GetGeometry()->GetLArTPCMap());
570 
571  if (larTPCMap.empty())
572  {
573  std::cout << "LArGeometryHelper::GetSigmaUVW - LArTPC description not registered with Pandora as required " << std::endl;
574  throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
575  }
576 
577  const LArTPC *const pFirstLArTPC(larTPCMap.begin()->second);
578  const float sigmaUVW(pFirstLArTPC->GetSigmaUVW());
579 
580  for (const LArTPCMap::value_type &mapEntry : larTPCMap)
581  {
582  const LArTPC *const pLArTPC(mapEntry.second);
583 
584  if (std::fabs(sigmaUVW - pLArTPC->GetSigmaUVW()) > maxSigmaDiscrepancy)
585  {
586  std::cout << "LArGeometryHelper::GetSigmaUVW - Plugin does not support provided LArTPC configurations " << std::endl;
587  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
588  }
589  }
590 
591  return sigmaUVW;
592 }
593 
594 } // namespace lar_content
pandora::CartesianVector GetGlobalMinLayerDirection() const
Get global direction corresponding to the fit result in minimum fit layer.
std::set< unsigned int > UIntSet
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
Header file for the cluster helper class.
Header file for the lar two dimensional sliding fit result class.
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
Double_t Z2
Definition: plot.C:263
Double_t Z1
Definition: plot.C:263
const pandora::Cluster * GetCluster() const
Get the address of the cluster, if originally provided.
pandora::StatusCode GetGlobalFitPositionAtX(const float x, pandora::CartesianVector &position) const
Get global fit position for a given input x coordinate.
HitType
Definition: HitType.h:12
pandora::CartesianVector GetGlobalMaxLayerDirection() const
Get global direction corresponding to the fit result in maximum fit layer.
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.