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 
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  {
25  Add(c);
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  //......................................................................
142  void Chain::Add(const Cand& c)
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)
153  Add(c);
154  else
155  AddFront(c);
156  }
157 
158  //......................................................................
159  void Chain::AddFront(const Cand& c)
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 
190  while(!chunks.empty() && chunks.back().AllDead()) chunks.pop_back();
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  }
235  hmap.Add(hit);
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)
271  if(!chunks[n].AllDead()) ++ret;
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()
Definition: readHits.C:15
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
void AddFront(const Cand &c)
Definition: Chain.cxx:159
Direction
Definition: Types.h:5
void Add(const Cand &c)
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: event.h:1
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