Chain.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file Chain.cxx
3 // \version \$Id: Chain.cxx,v 1.27 2012-11-01 06:05:03 bckhouse Exp \$
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6
8
9 #include "DiscreteTracker/Cand.h"
10
11
12
13
14 #include "RecoBase/CellHit.h"
15 #include "RecoBase/HitMap.h"
16 #include "RecoBase/Track.h"
17
18 #include "Utilities/func/MathUtil.h"
19
20 namespace dt
21 {
22  //......................................................................
24  {
26  }
27
28  //......................................................................
30  {
32
33  // Special case for vertical cands (who'll otherwise estimate themselves
34  // arbitrarily at a slope depending if they're marked "up" or "down").
35  if(ExtentPlane() == 0 && AllHits().size() > 1){
36  const Chunk& chunk = FirstCand().FirstChunk();
37  double xyz0[3];
38  geom->Plane(chunk.Plane())->Cell(chunk.TopHitCell())->GetCenter(xyz0);
39  double xyz1[3];
40  geom->Plane(chunk.Plane())->Cell(chunk.BottomHitCell())->GetCenter(xyz1);
41
44  xyz0[view], xyz0[2],
45  xyz1[view]-xyz0[view], xyz1[2]-xyz0[2]);
46
47  tr.AppendTrajectoryPoint(xyz1[view], xyz1[2]);
48  return tr;
49  }
50
51  // The algorithm is to take the weighted average of the y-positions of all
52  // relevant candidates at each z-position. (Convention is that the non-z
53  // coordinate is called "y").
54
55  // Estimate the start and end points of all the component cands upfront
56  std::vector<TVector3> p0s, p1s;
57  for(std::list<Cand>::const_iterator candIt = fCands.begin();
58  candIt != fCands.end(); ++candIt){
59  double z1, y1, z2, y2;
60  candIt->EstimateStraightLine(z1, y1, z2, y2);
61  p0s.push_back(TVector3(0, y1, z1));
62  p1s.push_back(TVector3(0, y2, z2));
63  }
64
65  std::vector<TVector3> pts;
66
67  // TODO, figure out where -1e10 Z's come from
68  const double minZ = std::max(0., p0s[0].Z());
69  const double maxZ = std::min(p1s.back().Z(), geom->DetLength());
70
71  const double kZStep = 5; // cm
72  const int stepMax = int((maxZ-minZ)/kZStep);
73
74  for(int step = 0; step <= stepMax; ++step){
75  const double z = minZ+(double(step)/stepMax)*(maxZ-minZ);
76
77  TVector3 meanp;
78  double totweight = 0;
79
80  for(unsigned int pIdx = 0; pIdx < p0s.size(); ++pIdx){
81  const TVector3& p0 = p0s[pIdx];
82  const TVector3& p1 = p1s[pIdx];
83  const double dz = p1.Z()-p0.Z();
84
85  const TVector3 p = p0*((p1.Z()-z)/dz) + p1*((z-p0.Z())/dz);
86
87  // Give more weight to cands near their centres
88  double weight = (z-p0.Z())/dz;
89  if(weight > 0.5) weight = 1-weight;
90  weight += 1e-3; // Avoid zeros at very ends
91  if(weight < 0) continue;
92
93  meanp += p*weight;
94  totweight += weight;
95  }
96
97  // This can happen when two cands are connected because they intersect at
98  // a shallow angle inside a large dead region. There's probably a neater
99  // way to illustrate this in the trajectory, but it's rare, so this will
100  // do.
101  if(totweight == 0) continue; // No cands cover this range
102
103  meanp *= 1/totweight;
104
105  pts.push_back(meanp);
106  } // end for step
107
108  // TODO
109  // assert(pts.size() >= 2);
110  assert(!pts.empty());
111  if(pts.size() == 1){
112  mf::LogWarning("DiscreteTracker") << "Only one trajectory point. Fabricating a second one." << std::endl;
113  pts.push_back(pts[0]+TVector3(0, 0, 1));
114  }
115
116  const TVector3 dir = pts[1]-pts[0];
117
118  rb::Track tr(AllHits(),
119  fCands.begin()->View(),
120  pts[0].Y(), pts[0].Z(), dir.Y(), dir.Z());
121
122  for(unsigned int n = 1; n < pts.size(); ++n)
123  tr.AppendTrajectoryPoint(pts[n].Y(), pts[n].Z());
124
125  return tr;
126  }
127
128  //......................................................................
129  std::vector<rb::Track> Chain::ToDebugTracks() const
130  {
131  std::vector<rb::Track> ret;
132
133  for(std::list<Cand>::const_iterator candIt = fCands.begin();
134  candIt != fCands.end(); ++candIt){
135  ret.push_back(candIt->ToTrack());
136  }
137
138  return ret;
139  }
140
141  //......................................................................
143  {
144  assert(c.NChunks() > 0);
145  assert(c.NSegs() > 1);
146  fCands.push_back(c);
147  }
148
149  //......................................................................
150  void Chain::Add(const Cand& c, Direction dir)
151  {
152  if(dir == kDownstream)
154  else
156  }
157
158  //......................................................................
160  {
161  assert(c.NChunks() > 0);
162  assert(c.NSegs() > 1);
163  fCands.push_front(c);
164  }
165
166  //......................................................................
167  unsigned int Chain::FirstPlane() const
168  {
169  return FirstCand().FirstChunk().Plane();
170  }
171
172  //......................................................................
173  unsigned int Chain::LastPlane() const
174  {
175  return LastCand().LastChunk().Plane();
176  }
177
178  //......................................................................
179  unsigned int Chain::ExtremalPlane(Direction dir) const
180  {
181  return ExtremalCand(dir).ExtremalPlane(dir);
182  }
183
184  //......................................................................
186  {
187  std::vector<Chunk> chunks = AllChunks();
188  if(dir == kUpstream) std::reverse(chunks.begin(), chunks.end());
189
191
192  return chunks.back().Plane();
193  }
194
195  //......................................................................
196  const Chunk& Chain::FirstChunk() const
197  {
198  return FirstCand().FirstChunk();
199  }
200
201  //......................................................................
202  const Chunk& Chain::LastChunk() const
203  {
204  return LastCand().LastChunk();
205  }
206
207  //......................................................................
208  const Segment& Chain::FirstSeg() const
209  {
210  return FirstCand().FirstSeg();
211  }
212
213  //......................................................................
214  const Segment& Chain::LastSeg() const
215  {
216  return LastCand().LastSeg();
217  }
218
219  //......................................................................
221  {
223
224  // To eliminate duplicates
225  rb::HitMap hmap;
226
227  for(std::list<Cand>::const_iterator it = fCands.begin();
228  it != fCands.end(); ++it){
229  const art::PtrVector<rb::CellHit>& hits = it->AllHits();
230  for(unsigned int hitIdx = 0; hitIdx < hits.size(); ++hitIdx){
231  const art::Ptr<rb::CellHit> hit = hits[hitIdx];
232  if(!hmap.CellExists(hit->Plane(), hit->Cell())){
233  ret.push_back(hit);
234  }
236  }
237  } // end for it
238
239  // TODO
240  // assert(!ret.empty());
241  return ret;
242  }
243
244  //......................................................................
245  std::vector<Chunk> Chain::AllChunks() const
246  {
247  // TODO: this doesn't take into account the possibility that chunks on the
248  // same plane don't entirely correspond (due to dead channels?)
249  std::map<unsigned int, Chunk> retMap;
250
251  for(std::list<Cand>::const_iterator it = fCands.begin(); it != fCands.end(); ++it){
252  const std::list<Chunk> chunks = it->AllChunks();
253  for(std::list<Chunk>::const_iterator it = chunks.begin(); it != chunks.end(); ++it){
254  retMap[it->Plane()] = *it;
255  }
256  }
257
258  // Come out sorted in plane order
259  std::vector<Chunk> ret;
260  for(std::map<unsigned int, Chunk>::iterator it = retMap.begin(); it != retMap.end(); ++it)
261  ret.push_back(it->second);
262  return ret;
263  }
264
265  //......................................................................
266  int Chain::NHitPlanes() const
267  {
268  std::vector<Chunk> chunks = AllChunks();
269  int ret = 0;
270  for(unsigned int n = 0; n < chunks.size(); ++n)
272
273  return ret;
274  }
275
276  //......................................................................
278  {
279  const unsigned int origSize = AllHits().size();
280
281  assert(origSize > 0);
282
283  while(true){
284  const Cand cand = fCands.front();
285  fCands.pop_front();
286  if(fCands.empty() || AllHits().size() < origSize){
287  fCands.push_front(cand);
288  break;
289  }
290  }
291
292  while(true){
293  const Cand cand = fCands.back();
294  fCands.pop_back();
295  if(fCands.empty() || AllHits().size() < origSize){
296  fCands.push_back(cand);
297  break;
298  }
299  }
300
301  fCands.front().TrimFront();
302  fCands.back().TrimBack();
303  }
304
305  //......................................................................
306  void Chain::Truncate(unsigned int plane, Direction end)
307  {
308  if(end == kUpstream){
309  while(FirstPlane() < plane){
310  fCands.begin()->PopFirstChunk();
311  if(FirstCand().AllChunks().empty()) fCands.pop_front();
312  }
313  while(FirstSeg().plane < int(plane)){
314  if(fCands.begin()->NSegs() > 1)
315  fCands.begin()->PopFirstSeg();
316  else
317  fCands.pop_front();
318  }
319  }
320  else{
321  while(LastPlane() > plane){
322  fCands.rbegin()->PopLastChunk();
323  if(LastCand().AllChunks().empty()) fCands.pop_back();
324  }
325  while(LastSeg().plane > int(plane)){
326  if(fCands.rbegin()->NSegs() > 1)
327  fCands.rbegin()->PopLastSeg();
328  else
329  fCands.pop_back();
330  }
331  }
332
333  if(!AllHits().empty()) TrimEnds();
334  }
335
336  //......................................................................
338  {
340
341  const std::vector<Chunk> chunks = AllChunks();
342  for(unsigned int n = 0; n < chunks.size(); ++n)
343  ret[chunks[n].Plane()].push_back(chunks[n]);
344
345  return ret;
346  }
347
348  //......................................................................
349  bool compareByLength(const Chain& a, const Chain& b)
350  {
351  return a.ExtentPlane() < b.ExtentPlane();
352  }
353
354  //......................................................................
355  bool compareByStart(const Chain& a, const Chain& b)
356  {
357  return a.FirstPlane() < b.FirstPlane();
358  }
359
360  //......................................................................
361  bool compareByEnd(const Chain& a, const Chain& b)
362  {
363  return a.LastPlane() < b.LastPlane();
364  }
365
366 } // namespace
T max(const caf::Proxy< T > &a, T b)
bool compareByLength(const Chain &a, const Chain &b)
Definition: Chain.cxx:349
int NHitPlanes() const
Definition: Chain.cxx:266
int Plane() const
Definition: Chunk.h:23
set< int >::iterator it
Float_t y1[n_points_granero]
Definition: compare.C:5
const Segment & FirstSeg() const
Definition: Chain.cxx:208
const Var weight
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
bool compareByStart(const Chain &a, const Chain &b)
Definition: Chain.cxx:355
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 FirstPlane() const
Definition: Chain.cxx:167
Definition: Chain.py:1
unsigned int ExtremalPlane(Direction dir) const
Definition: Cand.h:68
double DetLength() const
std::map< int, std::vector< Chunk > > ChunkMap
Definition: View.h:27
geo::View_t View() const
Definition: Cand.h:40
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const Chunk & FirstChunk() const
Definition: Cand.cxx:353
const Segment & LastSeg() const
Definition: Cand.cxx:371
const PlaneGeo * Plane(unsigned int i) const
Chain()
Definition: Chain.h:32
unsigned int NChunks() const
Definition: Cand.h:41
Provides efficient lookup of CellHits by plane and cell number.
Definition: HitMap.h:22
unsigned int LastPlane() const
Definition: Chain.cxx:173
Float_t Z
Definition: plot.C:38
unsigned short Cell() const
Definition: CellHit.h:40
Definition: Cand.cxx:23
std::vector< Chunk > AllChunks() const
Definition: Chain.cxx:245
const Segment & LastSeg() const
Definition: Chain.cxx:214
rb::Track ToTrack() const
Definition: Chain.cxx:29
void hits()
Calculation and representation of a straight line passing through several "segment" windows...
Definition: Cand.h:25
int TopHitCell() const
Definition: Chunk.h:28
double dz[NP][NC]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
const double a
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
const Chunk & LastChunk() const
Definition: Chain.cxx:202
std::vector< rb::Track > ToDebugTracks() const
Definition: Chain.cxx:129
int BottomHitCell() const
Definition: Chunk.h:29
const Segment & FirstSeg() const
Definition: Cand.cxx:365
unsigned int ExtremalHitPlane(Direction dir) const
Definition: Chain.cxx:185
int NSegs() const
Definition: Cand.h:85
Definition: Chain.cxx:159
Direction
Definition: Types.h:5
Definition: Chain.cxx:142
Definition: View.py:1
z
Definition: test.py:28
dt::View::ChunkMap AsChunkMap() const
Definition: Chain.cxx:337
size_type size() const
Definition: PtrVector.h:308
const Chunk & FirstChunk() const
Definition: Chain.cxx:196
const Chunk & LastChunk() const
Definition: Cand.cxx:359
art::PtrVector< rb::CellHit > AllHits() const
Definition: Chain.cxx:220
unsigned int ExtremalPlane(Direction dir) const
Definition: Chain.cxx:179
void TrimEnds()
Remove Cands that add no actual hits. After building but before writing.
Definition: Chain.cxx:277
TDirectory * dir
Definition: macro.C:5
Definition: structs.h:12
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
unsigned int ExtentPlane() const
Definition: Chain.h:43
const Cand & FirstCand() const
Definition: Chain.h:45
bool CellExists(unsigned int planeIdx, unsigned int cellIdx) const
Does the map contain any cell at this position?
Definition: HitMap.cxx:81
void Truncate(unsigned int plane, Direction end)
Definition: Chain.cxx:306
Window the line must pass through from (z,y0)-(z,y1)
Definition: Segment.h:24
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
const Cand & ExtremalCand(Direction dir) const
Definition: Chain.h:47
const Cand & LastCand() const
Definition: Chain.h:46
std::list< Cand > fCands
Definition: Chain.h:80
bool compareByEnd(const Chain &a, const Chain &b)
Definition: Chain.cxx:361
void Add(const art::Ptr< rb::CellHit > &chit)
Definition: HitMap.cxx:37