ElasticArms.cxx
Go to the documentation of this file.
1 ///
2 /// \file ElasticArms.cxx
3 /// \brief The basic elastic arms energy function
4 /// \author messier@indiana.edu
5 /// \version $Id: ElasticArms.cxx,v 1.5 2012-05-31 01:08:49 messier Exp $
6 ///
9 #include <iostream>
10 #include <cmath>
11 #include <cstdlib>
12 #include <cassert>
13 using namespace earms;
14 
15 //
16 // These math utils are included in the NOvA software but they are so
17 // small its easier to just define them then link to all the libraries
18 // required to link to the existing versions.
19 //
20 static double sqr(double x) { return x*x; }
21 static double dsqrtoline(double x0, double y0,
22  double x1, double y1,
23  double x2, double y2)
24 {
25  assert(x1!=x2 || y1!=y2);
26  double dx = x2-x1;
27  double dy = y2-y1;
28  return sqr(dx*(y1-y0)-dy*(x1-x0)) / (dx*dx+dy*dy);
29 }
30 
31 //......................................................................
32 
34  unsigned int m,
35  double l,
36  double lv) :
37  fProngFit(1),
38  fT(200.0),
39  fBeta(1.0/fT),
40  fBetaUnchanged(false),
41  fLambda(l),
42  fLambdaV(lv),
43  fN(n),
44  fXorY(n),
45  fView(n),
46  fZ(n),
47  fWsqr(n),
48  fDhi(n),
49  fDaa(m),
50  fMia(n),
51  fMiaUnchanged(n),
52  fMiaRowUnchanged(n),
53  fMiaColUnchanged(m),
54  fExpBetaMia(n),
55  fVia(n),
56  fX0(0),
57  fY0(0),
58  fZ0(0),
59  fM(m),
60  fdXds(m),
61  fdYds(m),
62  fdZds(m)
63 {
64  unsigned int i;
65  for (i=0; i<n; ++i) {
66  fMia[i]. resize(m);
67  fMiaUnchanged[i].resize(m);
68  fExpBetaMia[i].resize(m);
69  fVia[i]. resize(m);
70  }
71 }
72 
73 //......................................................................
74 
75 void ElasticArms::SetHit(unsigned int i,
76  double z,
77  double xory,
78  double sigma,
79  int view)
80 {
81  assert(i<fN);
82  fZ[i] = z;
83  fXorY[i] = xory;
84  fWsqr[i] = 1.0/(sigma*sigma);
85  fView[i] = view;
86 }
87 
88 //......................................................................
89 
90 void ElasticArms::SetT(double t)
91 {
92  fT = t;
93  fBeta = 1.0/fT;
94  fBetaUnchanged = false;
96 }
97 
98 //......................................................................
99 
100 void ElasticArms::SetLambda(double lam)
101 {
102  fLambda = lam;
104 }
105 
106 //......................................................................
107 
109 {
110  unsigned int i; // hit index
111  unsigned int a; // track index
112 
113  for(i = 0; i < fN; ++i) fMiaRowUnchanged[i] = true;
114  for(a = 0; i < fM; ++a) fMiaColUnchanged[a] = true;
115 
116  //
117  // Loop for each track a of the set of M track templates
118  //
119  for (a=0; a<fM; ++a) {
120 
121  // Update track end point
122  double trkX1 = fX0 + fdXds[a];
123  double trkY1 = fY0 + fdYds[a];
124  double trkZ1 = fZ0 + fdZds[a];
125 
126  //
127  // Loop for each hit i of the set of N hits
128  //
129  double trkXorY0;
130  double trkXorY1;
131  for (i=0; i<fN; ++i) {
132  double dot;
133  //
134  // Set coordinates based on which view the hit is in
135  //
136  switch (fView[i]) {
137  case geo::kX:
138  trkXorY0 = fX0;
139  trkXorY1 = trkX1;
140  dot = (trkZ1-fZ0)*(fZ[i]-fZ0) + (trkXorY1-fX0)*(fXorY[i]-fX0);
141  break;
142  case geo::kY:
143  trkXorY0 = fY0;
144  trkXorY1 = trkY1;
145  dot = (trkZ1-fZ0)*(fZ[i]-fZ0) + (trkXorY1-fY0)*(fXorY[i]-fY0);
146  break;
147  default: abort();
148  }
149  //
150  // If we're doing a line fit or it the hit is in the forward
151  // hemisphere use the distance from the hit to the arm template
152  //
153  const double oldmia = fMia[i][a];
154  if (fProngFit==0 || dot>0.0) {
155  fMia[i][a] = fWsqr[i]*dsqrtoline(fZ[i], fXorY[i],
156  fZ0, trkXorY0,
157  trkZ1, trkXorY1);
158  }
159  else {
160  //
161  // In the backwards hemisphere take the distance to the
162  // vertex. Make it get worse rapidly for differences larger
163  // than one sigma
164  //
165  fMia[i][a] = fWsqr[i]*(sqr(fZ[i]-fZ0)+sqr(fXorY[i]-trkXorY0));
166  if (fMia[i][a]>1.0) fMia[i][a] *= fMia[i][a];
167  }
168  const bool unchanged = (oldmia == fMia[i][a]);
169  fMiaUnchanged[i][a] = unchanged;
170  if(!unchanged){
171  fMiaRowUnchanged[i] = false;
172  fMiaColUnchanged[a] = false;
173  }
174 
175  // We only need to update fExpBetaMia if either fBeta or fMia has changed,
176  // and since it's expensive to call exp(), only do it if needed.
177  if(!fBetaUnchanged || !fMiaUnchanged[i][a]){
178  // If the argument of exp() is less than -746 in double precision, the
179  // result is zero. Saves 25% of all ElasticArms time to avoid these calls.
180  const double exparg = -fBeta*fMia[i][a];
181  fExpBetaMia[i][a] = exparg < -746? 0: exp(exparg);
182  }
183  } // loop on i hits
184  } // loop on a tracks
185 }
186 
187 //......................................................................
188 
190 {
191  static const double TINY = 1e-16;
192  unsigned int i;
193  unsigned int a;
194  unsigned int b;
195 
196  double d1, d2, denom;
197 
198  d1 = fExpBetaLambda+TINY;
199  for (i=0; i<fN; ++i) {
200  if(fBetaUnchanged && fMiaRowUnchanged[i]) continue;
201 
202  d2 = TINY;
203  for (b=0; b<fM; ++b) d2 += fExpBetaMia[i][b];
204 
205  denom = d1+d2;
206  for (a=0; a<fM; ++a) {
207  fVia[i][a] = fExpBetaMia[i][a]/denom;
208  }
209  } // loop on i=0...N hits
210 }
211 
212 //......................................................................
213 
215 {
216  unsigned int i, a;
217  double dminx, dminy;
218  double mavex, mavey;
219  double nx, ny;
220 
221  for (a=0; a<fM; ++a) {
222  if(fMiaColUnchanged[a]) continue;
223  //
224  // Compute the average weight for hits on this track
225  //
226  nx = 0.0;
227  ny = 0.0;
228  mavex = 0.0;
229  mavey = 0.0;
230  for (i=0; i<fN; ++i) {
231  switch (fView[i]) {
232  case 0: mavex += fVia[i][a]; nx += 1.0; break;
233  case 1: mavey += fVia[i][a]; ny += 1.0; break;
234  default: abort();
235  }
236  }
237  mavex /= nx;
238  mavey /= ny;
239 
240  dminx = 99E99;
241  dminy = 99E99;
242  for (i=0; i<fN; ++i) {
243  switch (fView[i]) {
244  case 0:
245  if (fVia[i][a]>=mavex && fDhi[i]<dminx) dminx = fDhi[i];
246  break;
247  case 1:
248  if (fVia[i][a]>=mavey && fDhi[i]<dminy) dminy = fDhi[i];
249  break;
250  default: abort();
251  }
252  }
253  fDaa[a] = sqrt(dminx)+sqrt(dminy);
254  }
255 }
256 
257 //......................................................................
258 
259 double ElasticArms::E()
260 {
261  unsigned int i;
262  unsigned int a;
263  double sViaMia = 0.0; // First term in the energy
264  double sVia = 0.0; // The inside sum in the second term
265  double ssVia = 0.0; // The outside sum in the second term
266 
267  //
268  // Update all the pieces that go into the calculation
269  //
270  this->UpdateMia();
271  this->UpdateVia();
272  this->UpdateDaa();
273  fBetaUnchanged = true;
274 
275  //
276  // Loop for each hit i of N
277  //
278  for (i=0; i<fN; ++i) {
279  //
280  // Loop for each track a of M
281  //
282  sVia = 0.0;
283  for (a=0; a<fM; ++a) {
284  sViaMia += fVia[i][a]*fMia[i][a];
285  sVia += fVia[i][a];
286  }
287  ssVia += sqr(sVia-1.0);
288  }
289 
290  //
291  // Loop for each prong a of M
292  //
293  double sDaa = 0.0;
294  for (a=0; a<fM; ++a) sDaa += fDaa[a];
295  sDaa *= 2.0/fLambdaV;
296 
297  //
298  // Three terms:
299  //
300  // 1st accounts for distance from hits to prongs
301  // 2nd is a penalty term for hits labeled as noise
302  // 3rd is a penalty term for arms that are detached from the vertex
303  //
304  return sViaMia + fLambda*ssVia + sDaa;
305 }
306 
307 //......................................................................
308 
309 void ElasticArms::SetArm(unsigned int a, double theta, double phi)
310 {
311  assert(a<fM);
312  fdXds[a] = sin(theta)*cos(phi);
313  fdYds[a] = sin(theta)*sin(phi);
314  fdZds[a] = cos(theta);
315 }
316 
317 //......................................................................
318 
319 void ElasticArms::SetVertex(double x, double y, double z)
320 {
321  fX0 = x;
322  fY0 = y;
323  fZ0 = z;
324 
325  //
326  // Update the table of distances from the vertex to the hits
327  //
328  unsigned int i;
329  for (i=0; i<fN; ++i) {
330  switch (fView[i]) {
331  case 0:
332  fDhi[i] = sqr(fXorY[i]-fX0) + sqr(fZ[i]-fZ0);
333  break;
334  case 1:
335  fDhi[i] = sqr(fXorY[i]-fY0) + sqr(fZ[i]-fZ0);
336  break;
337  default: abort();
338  }
339  }
340 }
341 
342 //......................................................................
343 
344 void ElasticArms::GetVertex(double& x, double& y, double& z) const
345 {
346  x = fX0;
347  y = fY0;
348  z = fZ0;
349 }
350 
351 //......................................................................
352 
353 double ElasticArms::SumVi(unsigned int a) const
354 {
355  //
356  // Play some defense by checking the bounds of a against the number
357  // of arms in the object. a cannot be larger than (# arms)-1
358  //
359  if (a>=this->fM) {
360  std::cerr << __FILE__ << ":" << __LINE__ << " Bad arm number "
361  << a << ">=" << this->fM << std::endl;
362  return -99.0; // An obviously bogus result
363  }
364  unsigned int i;
365  double sum = 0.0;
366  for (i=0; i<fN; ++i) {
367  sum += fVia[i][a];
368  }
369  return sum;
370 }
371 
372 //......................................................................
373 
374 unsigned int ElasticArms::LeadingArm() const
375 {
376  unsigned int a;
377  unsigned int amax = 0;
378  double w;
379  double wmax = 0;
380  for (a=0; a<fM; ++a) {
381  w = this->SumVi(a);
382  if (w>wmax) {
383  wmax = w;
384  amax = a;
385  }
386  }
387  return amax;
388 }
389 
390 ////////////////////////////////////////////////////////////////////////
double fBeta
1/T in units of 1/sigma^2
Definition: ElasticArms.h:100
double fLambda
Noise penalty term in units of sigma^2.
Definition: ElasticArms.h:102
unsigned int LeadingArm() const
std::vector< double > fWsqr
Combined weight for hit 1/(sigmaH^2 + sigmaMS^2)
Definition: ElasticArms.h:112
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
void SetLambda(double lam)
T sqrt(T number)
Definition: d0nt_math.hpp:156
double fY0
Vertex y location (cm)
Definition: ElasticArms.h:127
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::vector< bool > fMiaColUnchanged
Is a column of fMia changed in UpdateMia?
Definition: ElasticArms.h:120
OStream cerr
Definition: OStream.cxx:7
std::vector< double > fDhi
From hit to vertex location^2 [cm^2].
Definition: ElasticArms.h:113
double fZ0
Vertex z location (cm)
Definition: ElasticArms.h:128
double fLambdaV
Distance scale for hits to appear on (cm)
Definition: ElasticArms.h:104
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Elastic Arms vertex finding algorithm.
unsigned int fM
Number of arms.
Definition: ElasticArms.h:129
bool fBetaUnchanged
Whether fBeta is the same since last call to UpdateMia.
Definition: ElasticArms.h:101
std::vector< double > fXorY
Transverse hit positions (cm)
Definition: ElasticArms.h:109
void SetHit(unsigned int i, double z, double xory, double sigma, int view)
Definition: ElasticArms.cxx:75
std::vector< double > fdZds
Track dz/ds.
Definition: ElasticArms.h:132
double dy[NP][NC]
double SumVi(unsigned int a) const
double dx[NP][NC]
const double a
static double sqr(double x)
Definition: ElasticArms.cxx:20
void resize(T &x, std::vector< int > dims)
Definition: resize.hpp:41
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
void SetT(double t)
Definition: ElasticArms.cxx:90
void GetVertex(double &x, double &y, double &z) const
unsigned int fN
Data to fit.
Definition: ElasticArms.h:108
std::vector< double > fDaa
Distance from vertex to prong [cm].
Definition: ElasticArms.h:114
z
Definition: test.py:28
double sigma(TH1F *hist, double percentile)
int fProngFit
1 - treat arms as prongs in fit, 0 - treat as lines
Definition: ElasticArms.h:98
std::vector< double > fZ
Z locations of hits (cm)
Definition: ElasticArms.h:111
double dot(const std::vector< double > &x, const std::vector< double > &y)
Definition: dot.hpp:10
double fX0
Model parameters.
Definition: ElasticArms.h:126
T sin(T number)
Definition: d0nt_math.hpp:132
void SetArm(unsigned int i, double theta, double phi)
std::vector< bool > fMiaRowUnchanged
Is a row of fMia changed in UpdateMia?
Definition: ElasticArms.h:119
std::vector< double > fdXds
Track dx/ds.
Definition: ElasticArms.h:130
matrix fVia
Strength of association hit i to track a.
Definition: ElasticArms.h:122
std::vector< std::vector< bool > > fMiaUnchanged
Is fMia changed in UpdateMia?
Definition: ElasticArms.h:118
const hit & b
Definition: hits.cxx:21
ElasticArms(unsigned int n, unsigned int m, double lambda, double lambdav)
Definition: ElasticArms.cxx:33
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
std::vector< int > fView
Which view is the hit in?
Definition: ElasticArms.h:110
matrix fExpBetaMia
exp(-fBeta*Mia) - compute once used many time
Definition: ElasticArms.h:121
double fExpBetaLambda
exp(-fBeta*fLambda) - computed once used many
Definition: ElasticArms.h:103
std::vector< double > fdYds
Track dy/ds.
Definition: ElasticArms.h:131
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
static double dsqrtoline(double x0, double y0, double x1, double y1, double x2, double y2)
Definition: ElasticArms.cxx:21
Float_t w
Definition: plot.C:20
double fT
A measure of the "temperature" in units of sigma^2.
Definition: ElasticArms.h:99
void SetVertex(double x, double y, double z)
matrix fMia
Distance^2 hit i to track a (sigma^2)
Definition: ElasticArms.h:117