HitList3D.cxx
Go to the documentation of this file.
2 #include <cstdlib>
3 #include <cmath>
4 #include <iostream>
5 #include <algorithm>
6 using namespace bpfit;
7 
9 HitList3D::HitList3D(unsigned int n) : std::vector<Hit3D>(n) { }
10 
11 //......................................................................
12 
14  double x,
15  double y,
16  double z,
17  double t,
18  double sig_x,
19  double sig_y,
20  double sig_z,
21  double sig_t,
22  double w)
23 {
24  this->push_back(Hit3D(v,x,y,z,t,sig_x,sig_y,sig_z,sig_t,w));
25 }
26 
27 //......................................................................
28 
29 void HitList3D::SetHit(unsigned int i,
30  geo::View_t v,
31  double x,
32  double y,
33  double z,
34  double t,
35  double sig_x,
36  double sig_y,
37  double sig_z,
38  double sig_t,
39  double w)
40 {
41  Hit3D& h = (*this)[i];
42  h.fView = v;
43  h.fX = x;
44  h.fY = y;
45  h.fZ = z;
46  h.fT = t;
47  h.fSigX = sig_x;
48  h.fSigY = sig_y;
49  h.fSigZ = sig_z;
50  h.fW = w;
51 }
52 
53 //......................................................................
54 ///
55 /// Supply a default interpolation routine. This makes a staight line
56 /// fit to the (up to) 4 closest hits from the other view.
57 ///
58 /// This function assumes the hits are sorted by z for optimum CPU
59 /// performance.
60 ///
62 {
63  unsigned int i, i0, i1, i2;
64  unsigned int npts = 6; // Minimum number of hits to use in interpolation
65  double dz = 100.0; // Range in z to search for opposite view hits
66  double dzmax = 800.0; // Largest z to search over
67 
68  for (i0=0; i0<hl.size(); ++i0) {
69  if (&h==&hl[i0]) break;
70  }
71 
72  unsigned int n = 0; // Number of hits used for interpolation
73  double Sw=0, Swx=0, Swy=0, Swxy=0, Swx2=0, wz;
74  while (dz<dzmax) {
75  for (i1=i0; i1!=0; --i1) {
76  if (fabs(h.fZ-hl[i1].fZ)>dz) { ++i1; break; }
77  }
78  for (i2=i0; i2<hl.size(); ++i2) {
79  if (fabs(h.fZ-hl[i2].fZ)>dz) break;
80  }
81  n = 0;
82  Sw = Swx = Swy = Swxy = Swx2 = 0;
83  for (i=i1; i<i2; ++i) {
84  //
85  // Skip hits from the same view
86  //
87  if (h.fView==hl[i].fView) continue;
88  if (fabs(h.fZ-hl[i].fZ)<dz) {
89  ++n;
90  Sw += hl[i].fW;
91  wz = hl[i].fZ * hl[i].fW;
92  Swx += wz;
93  Swx2 += wz*hl[i].fZ;
94  if (hl[i].fView==geo::kX) {
95  Swy += hl[i].fW*hl[i].fX;
96  Swxy += wz*hl[i].fX;
97  }
98  else if (hl[i].fView==geo::kY) {
99  Swy += hl[i].fW*hl[i].fY;
100  Swxy += wz*hl[i].fY;
101  }
102  else abort();
103  }
104  }
105  //
106  // Test doneness
107  //
108  if (n>=npts) break; // Have enough points to interpolate
109  //
110  // Not done, open up the window and keep searching
111  //
112  dz *= 2;
113  }
114  //
115  // Use what points we found to make a linear fit. In our coordinate
116  // system "z" plays the role of x and "xy" plays the tole of y for a
117  // fit y=m*x+b
118  //
119  double d = Sw*Swx2-Swx*Swx;
120  if (d==0.0) {
121  //
122  // If d is zero we have no reliable information from the other
123  // view. Put the hit at the detector center with large error.
124  //
125  if (h.fView==geo::kX) { h.fY = 0.0; h.fSigY = 1000; }
126  else if (h.fView==geo::kY) { h.fX = 0.0; h.fSigX = 1000; }
127  return;
128  }
129  double m = (Sw*Swxy-Swx*Swy)/d;
130  double b = -(Swx*Swxy-Swx2*Swy)/d;
131 
132  if (h.fView==geo::kX) {
133  h.fY = m*h.fZ+b;
134  h.fSigY = sqrt(m*m*h.fSigZ*h.fSigZ + h.fSigX*h.fSigX);
135  }
136  else if (h.fView==geo::kY) {
137  h.fX = m*h.fZ+b;
138  h.fSigX = sqrt(m*m*h.fSigZ*h.fSigZ + h.fSigY*h.fSigY);
139  }
140  else abort();
141 }
142 
143 //......................................................................
144 
146 {
147  if (interp==0) {
148  this->SortByZ();
149  interp = gsDefaultInterp;
150  }
151 
152  // ...for the sake of readability
153  HitList3D& hitlist = (*this);
154 
155  unsigned int i;
156  for (i=0; i<this->size(); ++i) {
157  (*interp)(hitlist[i], hitlist);
158  }
159 }
160 
161 //......................................................................
162 
163 static bool gsSortHit3DByZ(const Hit3D& h1, const Hit3D& h2)
164 {
165  return h1.fZ<h2.fZ;
166 }
167 
168 //......................................................................
169 
170 static bool gsSortHit3DByY(const Hit3D& h1, const Hit3D& h2)
171 {
172  return h1.fY>h2.fY;
173 }
174 
175 //......................................................................
176 
178 {
179  sort(this->begin(), this->end(), gsSortHit3DByZ);
180 }
181 
182 //......................................................................
183 
185 {
186  sort(this->begin(), this->end(), gsSortHit3DByY);
187 }
188 
189 ////////////////////////////////////////////////////////////////////////
static bool gsSortHit3DByZ(const Hit3D &h1, const Hit3D &h2)
Definition: HitList3D.cxx:163
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
T sqrt(T number)
Definition: d0nt_math.hpp:156
double fSigX
Error on x.
Definition: HitList3D.h:39
Vertical planes which measure X.
Definition: PlaneGeo.h:28
void(* Interp3D)(Hit3D &h, HitList3D &hl)
Definition: HitList3D.h:55
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
double fW
Weight of hit.
Definition: HitList3D.h:43
void CompleteOrtho(Interp3D interp=0)
Definition: HitList3D.cxx:145
double dz[NP][NC]
base_types push_back(int_type())
Contain 3D hit information by extrapolating ortho view coordinate from surrounding points in the orth...
double fX
X coordinate of hit.
Definition: HitList3D.h:35
void gsDefaultInterp(Hit3D &h, HitList3D &hl)
Definition: HitList3D.cxx:61
Float_t d
Definition: plot.C:236
geo::View_t fView
View hit belongs in (x or y)
Definition: HitList3D.h:34
z
Definition: test.py:28
double fT
time coordinate of hit
Definition: HitList3D.h:38
TH1F * h2
Definition: plot.C:45
static bool gsSortHit3DByY(const Hit3D &h1, const Hit3D &h2)
Definition: HitList3D.cxx:170
TH1F * h1
void AddHit(geo::View_t, double x, double y, double z, double t, double sig_x, double sig_y, double sig_z, double sig_t, double w=1.0)
Definition: HitList3D.cxx:13
double fSigY
Error on y.
Definition: HitList3D.h:40
const hit & b
Definition: hits.cxx:21
double fY
Y coordinate of hit.
Definition: HitList3D.h:36
Essential 3D information for tracking.
Definition: HitList3D.h:58
Float_t w
Definition: plot.C:20
void SetHit(unsigned int i, geo::View_t, double x, double y, double z, double t, double sig_x, double sig_y, double sig_z, double sig_t, double w=1.0)
Definition: HitList3D.cxx:29
double fZ
Z coordinate of hit.
Definition: HitList3D.h:37
double fSigZ
Error on z.
Definition: HitList3D.h:41