View.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file View.cxx
3 // \version $Id: View.cxx,v 1.22 2012-10-31 00:05:29 bckhouse Exp $
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "DiscreteTracker/View.h"
8 
10 #include "Geometry/Geometry.h"
11 #include "RecoBase/CellHit.h"
12 
14 
15 namespace dt
16 {
17 
18  //......................................................................
19  View::View(art::PtrVector<rb::CellHit> chits, bool considerBadChans)
20  {
21  assert(!chits.empty());
22 
23  fGeoView = chits[0]->View();
24 
27 
28  fChans.resize(geom->NPlanes());
29 
30  fPlanes = geom->GetPlanesByView(fGeoView);
31 
32  rb::Cluster clust(chits); // Just to get the extents
33  // Clumsy technique to remove all the planes outside the hit range
34  const unsigned int minPlane = clust.MinPlane();
35  for(unsigned int p = 0; p < minPlane; ++p) fPlanes.erase(p);
36  const unsigned int maxPlane = clust.MaxPlane();
37  for(unsigned int p = maxPlane+1; p < geom->NPlanes(); ++p)
38  fPlanes.erase(p);
39 
40  // Have the box we save be one bigger than the extent of the hit cells (but
41  // not larger than the detector).
42  const unsigned int minCell = std::max(clust.MinCell(fGeoView), 1u)-1;
43  fMinCell = minCell;
44  const unsigned int maxCell = std::min(clust.MaxCell(fGeoView)+1, geom->Plane(minPlane)->Ncells()-1);
45 
46  for(std::set<unsigned int>::iterator it = fPlanes.begin();
47  it != fPlanes.end(); ++it){
48  const unsigned int plane = *it;
49 
50  // Are you running with a geometry inconsistent with what you generated
51  // with?
52  assert(plane < geom->NPlanes());
53 
54  fChans[plane].reserve(maxCell-minCell+1);
55 
56  for(unsigned int cell = minCell; cell <= maxCell; ++cell){
57  dt::Channel chan(plane, cell);
58  if(considerBadChans &&
59  (bcl->IsBad(plane, cell) || bcl->IsLowEfficiency(plane, cell))){
60  chan.SetDead();
61  }
62 
63  fChans[plane].push_back(chan);
64  }
65  } // end for it (plane)
66 
67  for(unsigned int n = 0; n < chits.size(); ++n){
68  const art::Ptr<rb::CellHit> chit = chits[n];
69  const int plane = chit->Plane();
70  const int cell = chit->Cell();
71  dt::Channel& chan = fChans[plane][cell-minCell];
72  // Behave correctly in the face of eg weird CalHit configuration, and of
73  // channels that have bad efficiency but somehow worked this one time.
74  if(chan.Type() != Channel::kDead ||
75  (considerBadChans && bcl->IsLowEfficiency(plane, cell))){
76  chan.SetHit(chit);
77  }
78  }
79 
80  // Don't consider any planes that are completely dead
81  std::vector<unsigned int> deadPlanes;
82  for(std::set<unsigned int>::iterator it = fPlanes.begin();
83  it != fPlanes.end(); ++it){
84  const unsigned int plane = *it;
85 
86  bool alldead = true;
87  for(unsigned int cell = 0; cell < fChans[plane].size(); ++cell){
88  if(fChans[plane][cell].Type() != Channel::kDead){
89  alldead = false;
90  break;
91  }
92  } // end for cell
93  if(alldead) deadPlanes.push_back(plane);
94  } // end for it (plane)
95 
96  for(unsigned int n = 0; n < deadPlanes.size(); ++n)
97  fPlanes.erase(deadPlanes[n]);
98  }
99 
100  //......................................................................
101  int View::PrevPlane(unsigned int plane) const
102  {
103  std::set<unsigned int>::iterator it = fPlanes.find(plane);
104  if(it == fPlanes.end()) return -1;
105  if(it == fPlanes.begin()) return -1;
106  --it;
107  return *it;
108  }
109 
110  //......................................................................
111  int View::NextPlane(unsigned int plane) const
112  {
113  std::set<unsigned int>::iterator it = fPlanes.find(plane);
114  if(it == fPlanes.end()) return -1;
115  ++it;
116  if(it == fPlanes.end()) return -1;
117  return *it;
118  }
119 
120  //......................................................................
121  int View::AdjacentPlane(unsigned int plane, Direction dir) const
122  {
123  if(dir == kUpstream) return PrevPlane(plane);
124  return NextPlane(plane);
125  }
126 
127  //......................................................................
129  {
130  if(Outside(plane, cell)){
131  dt::Channel outside(plane, cell);
132  outside.SetDead();
133  return outside;
134  }
135 
136  return fChans[plane][cell-fMinCell];
137  }
138 
139  //......................................................................
141  {
142  assert(!Outside(plane, cell));
143 
144  return fChans[plane][cell-fMinCell];
145  }
146 
147  //......................................................................
149  {
150  ChunkMap chunks_map;
151 
152  for(std::set<unsigned int>::iterator it = fPlanes.begin();
153  it != fPlanes.end(); ++it){
154  const int plane = *it;
155 
156  bool inside = false;
157  int begin = -1;
158 
159  const unsigned int cellMax = fMinCell+fChans[plane].size();
160  for(unsigned int cell = fMinCell; cell < cellMax; ++cell){
161  const Channel::Type_t type = Channel(plane, cell).Type();
162  if(!inside && (type == Channel::kHit || type == Channel::kDead)){
163  inside = true;
164  begin = cell;
165  }
166  if(inside && (type != Channel::kHit && type != Channel::kDead)){
167  inside = false;
168  std::vector<Chunk> chunks = MakeChunkCombos(plane, begin, cell);
169  chunks_map[plane].insert(chunks_map[plane].begin(),
170  chunks.begin(), chunks.end());
171  }
172  } // end for cell
173  if(inside){
174  std::vector<Chunk> chunks = MakeChunkCombos(plane, begin, cellMax);
175  chunks_map[plane].insert(chunks_map[plane].begin(),
176  chunks.begin(), chunks.end());
177  }
178  } // end for it (plane)
179 
180  return chunks_map;
181  }
182 
183  //......................................................................
184  std::vector<Chunk> View::MakeChunkCombos(int plane,
185  int begin, int end) const
186  {
187  std::vector<Chunk> ret;
188 
189  assert(begin < end);
190 
191  // Pairs of hit cell and extremal neighbouring dead cell, ordered low to high
192  std::vector<std::pair<int, int> > lows, highs;
193 
194  // Every hit cell that has a boundary below, and every one with a boundary above
195  for(int k = begin; k < end; ++k){
196  const Channel::Type_t here = Channel(plane, k).Type();
197  const Channel::Type_t below = Channel(plane, k-1).Type();
198  const Channel::Type_t above = Channel(plane, k+1).Type();
199 
200  if(here == Channel::kHit && (k == begin || below == Channel::kDead)){
201  int lower = k;
202  // Run the dead cell along until we find the end of the dead region
203  while(Channel(plane, lower-1).Type() == Channel::kDead && lower-1 >= begin) --lower;
204  lows.push_back(std::make_pair(lower, k));
205  }
206 
207  if(here == Channel::kHit && (k == end-1 || above == Channel::kDead)){
208  int higher = k;
209  while(Channel(plane, higher+1).Type() == Channel::kDead && higher+1 < end) ++higher;
210  highs.push_back(std::make_pair(k, higher));
211  }
212  } // end for k
213 
214 
215  // Every combination of those starting and stopping points
216  for(unsigned int i = 0; i < lows.size(); ++i){
217  for(unsigned int j = 0; j < highs.size(); ++j){
218  if(lows[i].second > highs[j].first) continue;
219 
221  for(int k = lows[i].second; k <= highs[j].first; ++k){
222  if(Channel(plane, k).Type() == Channel::kHit)
223  hits.push_back(Channel(plane, k).GetHit());
224  }
225 
226  for(bool up: {false, true}){
227  Chunk chunk;
228  if(up) chunk = Chunk(plane, lows[i].first, lows[i].second, highs[j].first, highs[j].second, up, hits);
229  else chunk = Chunk(plane, highs[j].first, highs[j].second, lows[i].first, lows[i].second, up, hits);
230 
231  ret.push_back(chunk);
232  } // end for up
233  } // end for i
234  } // end for j
235 
236 
237  // Also, make a chunk for every contiguous block of dead channels. A track
238  // can go through these without needing to hit the live cells
239  // adjacent. Hopefully the scoring system weighting real hits higher means
240  // this won't cause problems.
241  for(int i = begin; i < end; ++i){
242  if(Channel(plane, i).Type() != Channel::kDead) continue;
243  const int from = i;
244  for(; i < end; ++i){
245  if(Channel(plane, i).Type() != Channel::kDead) break;
246  }
247  const int to = i-1;
248 
249  for(bool up: {false, true}){
250  const Chunk chunk(plane, from, to, from, to, up,
252  ret.push_back(chunk);
253  }
254  }
255 
256  assert(!ret.empty());
257  return ret;
258  }
259 
260  //......................................................................
261  bool View::Outside(int plane, int cell) const
262  {
263  return (plane < 0 || cell < int(fMinCell) ||
264  plane >= int(fChans.size()) ||
265  cell-int(fMinCell) >= int(fChans[plane].size()));
266  }
267 
268 } // namespace
bool Outside(int plane, int cell) const
Definition: View.cxx:261
T max(const caf::Proxy< T > &a, T b)
void SetDead()
Definition: Channel.cxx:19
set< int >::iterator it
unsigned short Plane() const
Definition: CellHit.h:39
Sequence of contiguous hits and dead cells all on the same plane.
Definition: Chunk.h:17
const char * p
Definition: xmltok.h:285
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
dt::Channel Channel(int plane, int cell) const
Definition: View.cxx:128
Type_t Type() const
Definition: Channel.cxx:45
A collection of associated CellHits.
Definition: Cluster.h:47
std::vector< Chunk > MakeChunkCombos(int plane, int begin, int end) const
Definition: View.cxx:184
std::map< int, std::vector< Chunk > > ChunkMap
Definition: View.h:27
unsigned int MaxCell(geo::View_t view) const
Definition: Cluster.cxx:518
std::set< unsigned int > fPlanes
Definition: View.h:49
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const PlaneGeo * Plane(unsigned int i) const
int AdjacentPlane(unsigned int plane, Direction dir) const
Definition: View.cxx:121
unsigned int fMinCell
Definition: View.h:45
::xsd::cxx::tree::type type
Definition: Database.h:110
int NextPlane(unsigned int plane) const
Definition: View.cxx:111
Representation of the state of a single detector channel.
Definition: Channel.h:19
unsigned short Cell() const
Definition: CellHit.h:40
Definition: Cand.cxx:23
ChunkMap MakeChunks() const
Definition: View.cxx:148
void hits()
Definition: readHits.C:15
geo::View_t fGeoView
Definition: View.h:43
void SetHit(art::Ptr< rb::CellHit > chit)
Definition: Channel.cxx:25
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
View(art::PtrVector< rb::CellHit > chits, bool considerBadChans)
Definition: View.cxx:19
const double j
Definition: BetheBloch.cxx:29
Direction
Definition: Types.h:5
bool empty() const
Definition: PtrVector.h:336
size_type size() const
Definition: PtrVector.h:308
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
int PrevPlane(unsigned int plane) const
Definition: View.cxx:101
TDirectory * dir
Definition: macro.C:5
const int cellMax
unsigned int MinCell(geo::View_t view) const
Definition: Cluster.cxx:472
const std::set< unsigned int > & GetPlanesByView(View_t v=kXorY) const
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
unsigned int NPlanes() const
T min(const caf::Proxy< T > &a, T b)
bool IsBad(int plane, int cell)
std::vector< std::vector< dt::Channel > > fChans
Should access this via Channel function.
Definition: View.h:47
Encapsulate the geometry of one entire detector (near, far, ndos)
bool IsLowEfficiency(int plane, int cell)