Util.cxx
Go to the documentation of this file.
1 #include "LEM/Util.h"
2 
4 
5 #include "TVector3.h"
6 
7 #include "Geometry/Geometry.h"
9 #include "RecoBase/Cluster.h"
10 
11 namespace lem
12 {
13  //......................................................................
14  int MeanCell(geo::View_t view, const rb::Cluster& slice)
15  {
16  return (slice.MinCell(view)+slice.MaxCell(view))/2;
17  }
18 
19  //......................................................................
20  void DefaultVertex(const rb::Cluster& slice,
21  int& plane, int& cell, int& cellOtherView,
22  bool reverse)
23  {
25 
26  // Fallback when everything goes wrong.
27 
28  // Just put the vertex at the top of the slice in the middle
29 
30  plane = reverse ? slice.MaxPlane() : slice.MinPlane();
31  if(plane < 0 || plane >= int(geom->NPlanes())) plane = 0; // arghhh
32  const geo::View_t view = geom->Plane(plane)->View();
33  cell = MeanCell(view, slice);
34  cellOtherView = MeanCell(geo::View_t(1-view), slice);
35  }
36 
37  //......................................................................
38  geo::CellUniqueId FindCellRobust(double x, double y, double z)
39  {
41 
42  for(int trial = 0; trial < 4; ++trial){
43  // Figure out what cell the vertex was in. Step 3cm in z in case of
44  // failure. If that didn't work, it must be in plastic between cells,
45  // move laterally in x or y.
46  double offsetX = 0, offsetY = 0;
47  if(trial == 1) offsetX = -2;
48  if(trial == 2) offsetY = -2;
49  if(trial == 3) offsetX = offsetY = -2;
50 
51  const geo::CellUniqueId id = geom->CellId(x+offsetX, y+offsetY, z,
52  0, 0, 1, 3);
53  // We found one that works
54  if(id) return id;
55  }
56  return 0;
57  }
58 
59  //......................................................................
60  void VertexToPlaneAndCell(const TVector3 vtx, const rb::Cluster& slice,
61  int& plane, int& cell, int& cellOtherView,
62  bool reverse)
63  {
64  // Picked up by failure cases below
65  DefaultVertex(slice, plane, cell, cellOtherView, reverse);
66 
68 
69  const geo::CellUniqueId id = FindCellRobust(vtx.X(),
70  vtx.Y(),
71  vtx.Z());
72 
73  // Take defaults
74  if(!id) return;
75 
76  geom->IdToCell(id, &plane, &cell);
77 
78  const int nextPlane = geom->NextPlaneOtherView(plane);
79  // Take default in other view
80 
81  if(nextPlane == geo::kPLANE_NOT_FOUND) return;
82 
83  double xyz[3];
84  int ncells = geom->Plane(nextPlane)->Ncells();
85  if(cellOtherView >= ncells)
86  cellOtherView = ncells-1;
87 
88  geom->Plane(nextPlane)->Cell(cellOtherView)->GetCenter(xyz);
89 
90  // z position of plane in other view
91  double zOther = xyz[2];
92 
93  // Sometimes this goes wrong and we have to iterate a bit
94  while(true){
95  const geo::CellUniqueId idOther = FindCellRobust(vtx.X(),
96  vtx.Y(),
97  zOther);
98  // Take default in other view
99  if(!idOther) return;
100 
101  // The cell we found must be in the plane we wanted
102  int checkPlane;
103  geom->IdToCell(idOther, &checkPlane, &cellOtherView);
104  // If so, we're good and return the answer
105  if(checkPlane == nextPlane) return;
106  // Otherwise, walk 1cm at a time in the right direction
107  if(checkPlane < nextPlane) zOther += 1;
108  if(checkPlane > nextPlane) zOther -= 1;
109  }
110  }
111 } // namespace
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
void VertexToPlaneAndCell(const TVector3 vtx, const rb::Cluster &slice, int &plane, int &cell, int &cellOtherView, bool reverse)
Definition: Util.cxx:60
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
A collection of associated CellHits.
Definition: Cluster.h:47
unsigned int MaxCell(geo::View_t view) const
Definition: Cluster.cxx:518
const PlaneGeo * Plane(unsigned int i) const
PID
Definition: FillPIDs.h:14
void DefaultVertex(const rb::Cluster &slice, int &plane, int &cell, int &cellOtherView, bool reverse)
Definition: Util.cxx:20
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
int MeanCell(geo::View_t view, const rb::Cluster &slice)
Definition: Util.cxx:14
z
Definition: test.py:28
geo::CellUniqueId FindCellRobust(double x, double y, double z)
Steps around a bit. Can still return zero for total failure.
Definition: Util.cxx:38
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
unsigned int MinCell(geo::View_t view) const
Definition: Cluster.cxx:472
void geom(int which=0)
Definition: geom.C:163
int ncells
Definition: geom.C:124
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
unsigned int NPlanes() const
const CellGeo * IdToCell(const CellUniqueId &id, int *iplane, int *icell) const
Encapsulate the geometry of one entire detector (near, far, ndos)
const unsigned int NextPlaneOtherView(unsigned int p, int d=+1) const