Chunk.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file Chunk.cxx
3 // \version $Id: Chunk.cxx,v 1.15 2012-10-30 23:33:08 bckhouse Exp $
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
8 
9 #include "Geometry/Geometry.h"
11 #include "RecoBase/CellHit.h"
12 
14 
15 namespace dt
16 {
17  // From gdml file. Not exposed by the Geometry service
18  const double kCornerR = .94; // cm
19 
20  //......................................................................
21  Chunk::Chunk(int plane, int cellInLo, int cellInHi, int cellOutLo, int cellOutHi, bool up,
23  : fPlane(plane),
24  fEntryCellLo(cellInLo), fEntryCellHi(cellInHi),
25  fExitCellLo(cellOutLo), fExitCellHi(cellOutHi),
26  fUp(up), fShowerChunk(false),
27  fSegsDirty(true),
28  fHits(hits)
29  {
30  assert(cellInLo <= cellInHi);
31  assert(cellOutLo <= cellOutHi);
32  }
33 
34  //......................................................................
35  bool Chunk::Matches(const Chunk& ch) const
36  {
37  if(fPlane != ch.fPlane) return false;
38 
39  if(BottomExtremeCell() != ch.BottomExtremeCell()) return false;
40  if(TopExtremeCell() != ch.TopExtremeCell()) return false;
41 
42  if(BottomHitCell() != ch.BottomHitCell()) return false;
43  if(TopHitCell() != ch.TopHitCell()) return false;
44 
45  return true;
46  }
47 
48  //......................................................................
49  bool Chunk::Up() const
50  {
51  return fUp;
52  }
53 
54  //......................................................................
55  void Chunk::SetUp(bool up)
56  {
57  fSegsDirty = true;
58 
59  if(up != fUp){
62  }
63  fUp = up;
64  }
65 
66  //......................................................................
68  double icept) const
69  {
71 
72  const geo::PlaneGeo* geoplane = geom->Plane(fPlane);
73 
75 
76  for(unsigned int n = 0; n < fHits.size(); ++n){
77  const art::Ptr<rb::CellHit>& hit = fHits[n];
78 
79  double xyz[3];
80  const geo::CellGeo* geocell = geoplane->Cell(hit->Cell());
81  geocell->GetCenter(xyz);
82  const double z0 = xyz[2];
83  const double dz = geocell->HalfD();
84  const double z1 = z0-dz;
85  const double z2 = z0+dz;
86  const double v0 = xyz[hit->View()];
87  const double dv = geocell->HalfW();
88  const double vlo = v0-dv;
89  const double vhi = v0+dv;
90 
91  double v1 = icept+grad*z1;
92  double v2 = icept+grad*z2;
93  if(v1 > v2) std::swap(v1, v2);
94 
95  if(v1 <= vhi && v2 >= vlo) ret.push_back(hit);
96  } // end for n
97 
98  return ret;
99  }
100 
101  //......................................................................
102  std::pair<Segment, Segment> Chunk::GetSegs() const
103  {
104  if(!fSegsDirty) return fSegs;
105  fSegsDirty = false;
106 
108 
109  const geo::PlaneGeo* plane = geom->Plane(fPlane);
110 
111  double v0lo, v1lo, v2lo, v3lo, z1lo, z2lo;
112  CellPoints(plane, BottomExtremeCell(),
113  v0lo, v1lo, v2lo, v3lo,
114  z1lo, z2lo);
115 
116  double v0hi, v1hi, v2hi, v3hi, z1hi, z2hi;
117  CellPoints(plane, TopExtremeCell(),
118  v0hi, v1hi, v2hi, v3hi,
119  z1hi, z2hi);
120 
121  double junk;
122  if(BottomHitCell() != BottomExtremeCell()){
123  CellPoints(plane, BottomHitCell(),
124  junk, junk, v2lo, v3lo,
125  junk, junk);
126  }
127  if(TopHitCell() != TopExtremeCell()){
128  CellPoints(plane, TopHitCell(),
129  v0hi, v1hi, junk, junk,
130  junk, junk);
131  }
132 
133  // If this is a completely dead chunk, it's OK to enter in the crack part
134  // on one side and exit in that same part without ever traversing any real
135  // cell.
136  if(fHits.empty()){
137  v1lo = v0lo;
138  v1hi = v0hi;
139  v2lo = v3lo;
140  v2hi = v3hi;
141  }
142 
143  // If the track is upgoing it can enter the plane in plastic below the
144  // lowest cell, and exit in the plastic above the highest cell. But it
145  // mustn't enter higher than the scintillator of the first cell, or lower
146  // than the scintillator of the last cell. Vice versa for downgoing.
147 
148  if(fShowerChunk){
149  if(fUp){
150  fSegs.first = Segment(z1lo, v0lo-kCornerR, v2hi,
151  z1lo+kCornerR, v0lo, v3hi,
153  fSegs.second = Segment(z2hi-kCornerR, v0lo, v3hi,
154  z2hi, v1lo, v3hi+kCornerR,
155  false, fPlane, fExitCellLo, fExitCellHi);
156  }
157  else{
158  fSegs.first = Segment(z1hi, v1lo, v3hi+kCornerR,
159  z1hi+kCornerR, v0lo, v3hi,
161  fSegs.second = Segment(z2lo-kCornerR, v0lo, v3hi,
162  z2lo, v0lo-kCornerR, v2hi,
163  false, fPlane, fExitCellLo, fExitCellHi);
164  }
165 
166  return fSegs;
167  }
168 
169  if(fUp){
170  // bottom left
171  fSegs.first = Segment(z1lo, v0lo-kCornerR, v2lo,
172  z1lo+kCornerR, v0lo, v3hi,
174  // top right
175  fSegs.second = Segment(z2hi-kCornerR, v0lo, v3hi,
176  z2hi, v1hi, v3hi+kCornerR,
177  false, fPlane, fExitCellLo, fExitCellHi);
178  }
179  else{
180  // top left
181  fSegs.first = Segment(z1hi, v1hi, v3hi+kCornerR,
182  z1hi+kCornerR, v0lo, v3hi,
184  // bottom right
185  fSegs.second = Segment(z2lo-kCornerR, v0lo, v3hi,
186  z2lo, v0lo-kCornerR, v2lo,
187  false, fPlane, fExitCellLo, fExitCellHi);
188  }
189 
190  return fSegs;
191  }
192 
193  //......................................................................
195  {
196  // Single hit
197  if(fEntryCellLo == fExitCellHi &&
198  fEntryCellHi == fExitCellLo) return false;
199  return true;
200  }
201 
202  //......................................................................
204  double& v0, double& v1, double& v2, double& v3,
205  double& z0, double& z1) const
206  {
207  const geo::View_t view = plane->View();
208 
209  // Measurements of the cell itself
210  double vc, dv, zc, dz;
211  CellMeasurements(view, plane->Cell(cell), vc, dv, zc, dz);
212  v1 = vc-dv;
213  v2 = vc+dv;
214  z0 = zc-dz;
215  z1 = zc+dz;
216 
217  // Would handle bad channels by walking these until the channel they're on
218  // is good.
219  const geo::CellGeo* cellabove = plane->Cell(cell+1);
220  const geo::CellGeo* cellbelow = plane->Cell(cell-1);
221 
222  double junk;
223  if(cellbelow){
224  // Lowest point not in another cell
225  CellMeasurements(view, cellbelow, v0, dv, junk, junk);
226  v0 += dv;
227  }
228  else{
229  // No cell below, must be edge of the detector, could be anywhere low
230  v0 = -1e3;
231  }
232 
233  if(cellabove){
234  // Highest point not in another cell
235  CellMeasurements(view, cellabove, v3, dv, junk, junk);
236  v3 -= dv;
237  }
238  else{
239  // No cell above, must be edge of the detector, could be anywhere high
240  v3 = +1e3;
241  }
242  }
243 
244  //......................................................................
246  double& v0, double& dv,
247  double& z0, double& dz) const
248  {
249  double xyz[3];
250  cell->GetCenter(xyz);
251  v0 = xyz[view];
252  dv = cell->HalfW();
253  z0 = xyz[2];
254  dz = cell->HalfD();
255  }
256 
257 } // namespace
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
void CellPoints(const geo::PlaneGeo *plane, int cell, double &v0, double &v1, double &v2, double &v3, double &z0, double &z1) const
Definition: Chunk.cxx:203
bool fShowerChunk
"Shower chunk" can enter and exit anywhere
Definition: Chunk.h:63
std::pair< Segment, Segment > GetSegs() const
Definition: Chunk.cxx:102
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Sequence of contiguous hits and dead cells all on the same plane.
Definition: Chunk.h:17
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
double HalfW() const
Definition: CellGeo.cxx:191
int fEntryCellLo
Definition: Chunk.h:61
const PlaneGeo * Plane(unsigned int i) const
int BottomExtremeCell() const
Definition: Chunk.h:31
bool Up() const
Definition: Chunk.cxx:49
art::PtrVector< rb::CellHit > fHits
Definition: Chunk.h:68
int fEntryCellHi
Definition: Chunk.h:61
static void grad(vari *vi)
Definition: grad.hpp:30
bool fSegsDirty
Definition: Chunk.h:65
int TopExtremeCell() const
Definition: Chunk.h:30
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
unsigned short Cell() const
Definition: CellHit.h:40
Definition: Cand.cxx:23
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
void hits()
Definition: readHits.C:15
int TopHitCell() const
Definition: Chunk.h:28
bool Matches(const Chunk &ch) const
Definition: Chunk.cxx:35
double dz[NP][NC]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
bool fUp
Definition: Chunk.h:62
art::PtrVector< rb::CellHit > HitsOnLine(double grad, double icept) const
Definition: Chunk.cxx:67
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
int fPlane
Definition: Chunk.h:60
bool IsWorthShowering() const
Definition: Chunk.cxx:194
int BottomHitCell() const
Definition: Chunk.h:29
void CellMeasurements(geo::View_t view, const geo::CellGeo *cell, double &v0, double &dv, double &z0, double &dz) const
Definition: Chunk.cxx:245
bool empty() const
Definition: PtrVector.h:336
size_type size() const
Definition: PtrVector.h:308
void SetUp(bool up)
Definition: Chunk.cxx:55
Definition: structs.h:12
std::pair< Segment, Segment > fSegs
Definition: Chunk.h:66
void geom(int which=0)
Definition: geom.C:163
const double kCornerR
Definition: Chunk.cxx:18
assert(nhit_max >=nhit_nbins)
Window the line must pass through from (z,y0)-(z,y1)
Definition: Segment.h:24
Chunk()
Definition: Chunk.h:20
Encapsulate the cell geometry.
Definition: CellGeo.h:25
int fExitCellHi
Definition: Chunk.h:61
Encapsulate the geometry of one entire detector (near, far, ndos)
int fExitCellLo
Definition: Chunk.h:61