Belle II Software  release-05-01-25
Lpav.cc
1 //-----------------------------------------------------------------------------
2 // $Id$
3 //-----------------------------------------------------------------------------
4 // Filename : Lpav.cc
5 // Section : TRG CDC
6 // Owner : Yoshihito Iwasaki
7 // Email : yoshihito.iwasaki@kek.jp
8 //-----------------------------------------------------------------------------
9 // Description :
10 //-----------------------------------------------------------------------------
11 // $Log$
12 //-----------------------------------------------------------------------------
13 
14 #include <cmath>
15 #include <iostream>
16 
17 #include "CLHEP/Vector/Sqr.h"
18 #include "trg/cdc/Lpav.h"
19 
20 namespace Belle2 {
26 //
27 // constants, enums and typedefs
28 //
30  extern "C" {
31  float prob_(float*, int*);
32  }
33 
35  static double err_dis_inv(double x, double y, double w, double a, double b)
36  {
37  if (a == 0 && b == 0) {
38  return w;
39  } else {
40  double f = x * b - y * a;
41  double rsq = x * x + y * y;
42  f *= f;
43  return w * rsq / f;
44  }
45  }
46 
47 //
48 // static data member definitions
49 //
50 
51 //
52 // constructors and destructor
53 //
54  TRGCDCLpav::TRGCDCLpav() :
55  m_wsum(),
56  m_xsum(),
57  m_ysum(),
58  m_xxsum(),
59  m_yysum(),
60  m_xysum(),
61  m_xrrsum(),
62  m_yrrsum(),
63  m_rrrrsum(),
64  m_wsum_temp(),
65  m_xav(),
66  m_yav(),
67  m_xyavp(),
68  m_rscale(),
69  m_xxavp(),
70  m_yyavp(),
71  m_xrravp(),
72  m_yrravp(),
73  m_rrrravp(),
74  m_sinrot(),
75  m_cosrot(),
76  m_nc(),
77  m_chisq() // 2019/07/31 by ytlai
78  {
79  clear();
80  }
81 
82 // TRGCDCLpav::TRGCDCLpav( const TRGCDCLpav& )
83 // {
84 // }
85 
86  TRGCDCLpav::~TRGCDCLpav()
87  {
88  }
89 
90 //
91 // assignment operators
92 //
93 // const TRGCDCLpav& TRGCDCLpav::operator=( const TRGCDCLpav& )
94 // {
95 // }
96 
97 //
98 // comparison operators
99 //
100 // bool TRGCDCLpav::operator==( const TRGCDCLpav& ) const
101 // {
102 // }
103 
104 // bool TRGCDCLpav::operator!=( const TRGCDCLpav& ) const
105 // {
106 // }
107 
108 //
109 // member functions
110 //
111  void TRGCDCLpav::calculate_average(double xi, double yi, double wi)
112  {
113  if (m_wsum <= 0) return;
114  m_wsum_temp = m_wsum + wi;
115  double rri(xi * xi + yi * yi);
116  double wrri(wi * rri);
117  double wsum_inv(1 / m_wsum_temp);
118  m_xav = (m_xsum + wi * xi) * wsum_inv;
119  m_yav = (m_ysum + wi * yi) * wsum_inv;
120 
121  double xxav((m_xxsum + wi * xi * xi) * wsum_inv);
122  double yyav((m_yysum + wi * yi * yi) * wsum_inv);
123  double xyav((m_xysum + wi * xi * yi) * wsum_inv);
124  double xrrav((m_xrrsum + xi * wrri) * wsum_inv);
125  double yrrav((m_yrrsum + yi * wrri) * wsum_inv);
126  double rrrrav((m_rrrrsum + wrri * rri) * wsum_inv);
127 
128  calculate_average_n(xxav, yyav, xyav, xrrav, yrrav, rrrrav);
129 
130  }
131 
132  void TRGCDCLpav::calculate_average(void)
133  {
134  if (m_wsum <= 0) return;
135  m_wsum_temp = m_wsum;
136  double wsum_inv(1 / m_wsum_temp);
137  m_xav = m_xsum * wsum_inv;
138  m_yav = m_ysum * wsum_inv;
139 
140  double xxav(m_xxsum * wsum_inv);
141  double yyav(m_yysum * wsum_inv);
142  double xyav(m_xysum * wsum_inv);
143  double xrrav(m_xrrsum * wsum_inv);
144  double yrrav(m_yrrsum * wsum_inv);
145  double rrrrav(m_rrrrsum * wsum_inv);
146 
147  calculate_average_n(xxav, yyav, xyav, xrrav, yrrav, rrrrav);
148  }
149 
150  void TRGCDCLpav::calculate_average_n(double xxav, double yyav, double xyav,
151  double xrrav, double yrrav, double rrrrav)
152  {
153  double xxav_p = xxav - m_xav * m_xav;
154  double yyav_p = yyav - m_yav * m_yav;
155  double xyav_p = xyav - m_xav * m_yav;
156  double rrav_p = xxav_p + yyav_p;
157 
158  double a = std::fabs(xxav_p - yyav_p);
159  double b = 4 * xyav_p * xyav_p;
160  double asqpb = a * a + b;
161  double rasqpb = std::sqrt(asqpb);
162  double splus = 1 + a / rasqpb;
163  double sminus = b / (asqpb * splus);
164  splus = std::sqrt(0.5 * splus);
165  sminus = std::sqrt(0.5 * sminus);
166 //C
167 //C== First require : SIGN(C**2 - S**2) = SIGN(XXAV - YYAV)
168 //C
169  if (xxav_p <= yyav_p) {
170  m_cosrot = sminus;
171  m_sinrot = splus;
172  } else {
173  m_cosrot = splus;
174  m_sinrot = sminus;
175  }
176 //C
177 //C== Require : SIGN(S) = SIGN(XYAV)*SIGN(C) (Assuming SIGN(C) > 0)
178 //C
179  if (xyav_p < 0) m_sinrot = - m_sinrot;
180 //*
181 //* We now have the smallest angle that guarantees <X**2> > <Y**2>
182 //*
183 //* To get the SIGN of the charge right, the new X-AXIS must point
184 //* outward from the orgin. We are free to change signs of both
185 //* COSROT and SINROT simultaneously to accomplish this.
186 //*
187 //* Choose SIGN of C wisely to be able to get the sign of the charge
188 //*
189  if (m_cosrot * m_xav + m_sinrot * m_yav <= 0) {
190  m_cosrot = - m_cosrot;
191  m_sinrot = - m_sinrot;
192  }
193  m_rscale = std::sqrt(rrav_p);
194  double cos2 = m_cosrot * m_cosrot;
195  double sin2 = m_sinrot * m_sinrot;
196  double cs2 = 2 * m_sinrot * m_cosrot;
197  double rrav_p_inv(1 / rrav_p);
198  m_xxavp = (cos2 * xxav_p + cs2 * xyav_p + sin2 * yyav_p) * rrav_p_inv;
199  m_yyavp = (cos2 * yyav_p - cs2 * xyav_p + sin2 * xxav_p) * rrav_p_inv;
200 
201  double xav2 = m_xav * m_xav;
202  double yav2 = m_yav * m_yav;
203  double xrrav_p = (xrrav - 2 * xxav * m_xav + xav2 * m_xav -
204  2 * xyav * m_yav + m_xav * yav2) - m_xav * rrav_p;
205  double yrrav_p = (yrrav - 2 * yyav * m_yav + yav2 * m_yav -
206  2 * xyav * m_xav + m_yav * xav2) - m_yav * rrav_p;
207  m_xrravp = (m_cosrot * xrrav_p + m_sinrot * yrrav_p) * rrav_p_inv / m_rscale;
208  m_yrravp = (- m_sinrot * xrrav_p + m_cosrot * yrrav_p) * rrav_p_inv / m_rscale;
209 
210  double rrav = xxav + yyav;
211  double rrrrav_p = rrrrav
212  - 2 * m_yav * yrrav - 2 * m_xav * xrrav
213  + rrav * (xav2 + yav2)
214  - 2 * m_xav * xrrav_p - xav2 * rrav_p
215  - 2 * m_yav * yrrav_p - yav2 * rrav_p;
216  m_rrrravp = rrrrav_p * rrav_p_inv * rrav_p_inv;
217  m_xyavp = 0;
218  }
219 
220  void TRGCDCLpav::calculate_average3(double xi, double yi, double wi)
221  {
222  if (m_wsum <= 0) return;
223  m_wsum_temp = m_wsum + wi;
224  double wsum_inv(1 / m_wsum_temp);
225  double rri(xi * xi + yi * yi);
226  m_xav = (m_xsum + wi * xi) * wsum_inv;
227  m_yav = (m_ysum + wi * yi) * wsum_inv;
228 
229  m_rscale = 1;
230  m_cosrot = 1;
231  m_sinrot = 0;
232  m_xxavp = (m_xxsum + wi * xi * xi) * wsum_inv;
233  m_xyavp = (m_xysum + wi * xi * yi) * wsum_inv;
234  m_yyavp = (m_yysum + wi * yi * yi) * wsum_inv;
235  double wrri(wi * rri);
236  m_xrravp = (m_xrrsum + xi * wrri) * wsum_inv;
237  m_yrravp = (m_yrrsum + yi * wrri) * wsum_inv;
238  m_rrrravp = (m_rrrrsum + rri * wrri) * wsum_inv;
239  }
240 
241  void TRGCDCLpav::calculate_average3(void)
242  {
243  if (m_wsum <= 0) return;
244  m_wsum_temp = m_wsum;
245  double wsum_inv(1 / m_wsum_temp);
246  m_xav = m_xsum * wsum_inv;
247  m_yav = m_ysum * wsum_inv;
248 
249  m_rscale = 1;
250  m_cosrot = 1;
251  m_sinrot = 0;
252  m_xxavp = m_xxsum * wsum_inv;
253  m_xyavp = m_xysum * wsum_inv;
254  m_yyavp = m_yysum * wsum_inv;
255  m_xrravp = m_xrrsum * wsum_inv;
256  m_yrravp = m_yrrsum * wsum_inv;
257  m_rrrravp = m_rrrrsum * wsum_inv;
258  }
259 
260 
261 //
262 // const member functions
263 //
264 
265 //
266 // static member functions
267 //
268 
270  std::ostream& operator<<(std::ostream& o, const TRGCDCLpav& a)
271  {
272 // o << "wsum=" << a.m_wsum << " xsum=" << a.m_xsum << " ysum=" << a.m_ysum
273 // << " xxsum=" << a.m_xxsum << " xysum=" << a.m_xysum
274 // << " yysum=" << a.m_yysum
275 // << " xrrsum=" << a.m_xrrsum << " yrrsum=" << a.m_yrrsum
276 // << " rrrrsum=" << a.m_rrrrsum;
277 // o << " rscale=" << a.m_rscale
278 // << " xxavp=" << a.m_xxavp << " yyavp=" << a.m_yyavp
279 // << " xrravp=" << a.m_xrravp << " yrravp=" << a.m_yrravp
280 // << " rrrravp=" << a.m_rrrravp << " cosrot=" << a.m_cosrot
281 // << " sinrot=" << a.m_sinrot
282 // << std::endl;
283  o << " nc=" << a.m_nc << " chisq=" << a.m_chisq << " " << (TRGCDCLpar&) a;
284  return o;
285  }
286 
287  double TRGCDCLpav::solve_lambda(void)
288  {
289  if (m_rscale <= 0) return -1;
290  double xrrxrr = m_xrravp * m_xrravp;
291  double yrryrr = m_yrravp * m_yrravp;
292  double rrrrm1 = m_rrrravp - 1;
293  double xxyy = m_xxavp * m_yyavp;
294 
295  double c0 = rrrrm1 * xxyy - xrrxrr * m_yyavp - yrryrr * m_xxavp;
296  double c1 = - rrrrm1 + xrrxrr + yrryrr - 4 * xxyy;
297  double c2 = 4 + rrrrm1 - 4 * xxyy;
298  double c4 = - 4;
299 //
300 //C COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS
301 //
302  double c2d = 2 * c2;
303  double c4d = 4 * c4;
304 //
305  double lambda = 0;
306 
307  double chiscl = m_wsum_temp * m_rscale * m_rscale;
308  double dlamax = 0.001 / chiscl;
309  const int ntry = 5;
310  int itry = 0;
311  double dlambda = dlamax;
312  while (itry < ntry && std::fabs(dlambda) >= dlamax) {
313  double cpoly = c0 + lambda * (c1 + lambda *
314  (c2 + lambda * lambda * c4));
315  double dcpoly = c1 + lambda * (c2d + lambda * lambda * c4d);
316  dlambda = - cpoly / dcpoly;
317  lambda += dlambda;
318  itry ++;
319  }
320  lambda = lambda < 0 ? 0 : lambda;
321  return lambda;
322  }
323 
324  double TRGCDCLpav::solve_lambda3(void)
325  {
326  if (m_rscale <= 0) return -1;
327  double xrrxrr = m_xrravp * m_xrravp;
328  double yrryrr = m_yrravp * m_yrravp;
329 //cnv double rrrrm1 = m_rrrravp - 1;
330 //cnv double xxyy = m_xxavp * m_yyavp;
331 
332  double a = m_rrrravp;
333  double b = xrrxrr + yrryrr - m_rrrravp * (m_xxavp + m_yyavp);
334  double c = m_rrrravp * m_xxavp * m_yyavp
335  - m_yyavp * xrrxrr - m_xxavp * yrryrr
336  + 2 * m_xyavp * m_xrravp * m_yrravp - m_rrrravp * m_xyavp * m_xyavp;
337  if (c >= 0 && b <= 0) {
338  return (-b - std::sqrt(b * b - 4 * a * c)) / 2 / a;
339  } else if (c >= 0 && b > 0) {
340  std::cout << " returning " << -1 << std::endl;
341  return -1;
342  } else if (c < 0) {
343  return (-b + std::sqrt(b * b - 4 * a * c)) / 2 / a;
344  }
345  return -1;
346  }
347 
348  double TRGCDCLpav::calculate_lpar(void)
349  {
350  double lambda = solve_lambda();
351 // changed on Oct-13-93
352 // if (lambda<=0) return -1;
353  if (lambda < 0) return -1;
354  double h11 = m_xxavp - lambda;
355  double h22 = m_yyavp - lambda;
356  if (h11 == 0.0) return -1;
357  double h14 = m_xrravp;
358  double h24 = m_yrravp;
359  double h34 = 1 + 2 * lambda;
360  double rootsq = (h14 * h14 / h11 / h11) + 4 * h34;
361  if (std::fabs(h22) > std::fabs(h24)) {
362  if (h22 == 0.0) return -1;
363  double ratio = h24 / h22;
364  rootsq += ratio * ratio ;
365  m_kappa = 1 / std::sqrt(rootsq);
366  m_beta = - ratio * m_kappa;
367  } else {
368  if (h24 == 0.0) return -1;
369  double ratio = h22 / h24;
370  rootsq = 1 + ratio * ratio * rootsq;
371  m_beta = 1 / std::sqrt(rootsq);
372  m_beta = h24 > 0 ? -m_beta : m_beta;
373  m_kappa = -ratio * m_beta;
374  }
375  m_alpha = - (h14 / h11) * m_kappa;
376  m_gamma = - h34 * m_kappa;
377 // if (lambda<0.0001) {
378 // std::cout << " lambda=" << lambda << " h34=" << h34
379 // << " rootsq=" << rootsq << " h22=" << h22
380 // << " h11=" << h11 << " h14=" << h14 << " h24=" << h24 <<
381 // " " << *this << std::endl;
382 // }
383 //
384 //C TRANSFORM THESE INTO THE LAB COORDINATE SYSTEM
385 //
386 //C FIRST GET KAPPA AND GAMMA BACK TO REAL DIMENSIONS
387 //
388  scale(m_rscale);
389 //
390 //C NEXT ROTATE ALPHA AND BETA
391 //
392  rotate(m_cosrot, -m_sinrot);
393 //
394 //C THEN TRANSLATE BY (XAV,YAV)
395 //
396  move(-m_xav, -m_yav);
397  if (m_yrravp < 0) neg();
398  if (lambda >= 0) m_chisq = lambda * m_wsum_temp * m_rscale * m_rscale;
399  return lambda;
400  }
401 
402  double TRGCDCLpav::calculate_lpar3(void)
403  {
404  double lambda = solve_lambda3();
405 // changed on Oct-13-93
406 // if (lambda<=0) return -1;
407  if (lambda < 0) return -1;
408  double h11 = m_xxavp - lambda;
409  double h22 = m_yyavp - lambda;
410  double h14 = m_xrravp;
411  double h24 = m_yrravp;
412  m_gamma = 0;
413  double h12 = m_xyavp;
414  double det = h11 * h22 - h12 * h12;
415  if (det != 0) {
416  double r1 = (h14 * h22 - h24 * h12) / (det);
417  double r2 = (h24 * h11 - h14 * h12) / (det);
418  double kinvsq = r1 * r1 + r2 * r2;
419  m_kappa = std::sqrt(1 / kinvsq);
420  if (h11 != 0) m_alpha = -m_kappa * r1;
421  else m_alpha = 1;
422  if (h22 != 0) m_beta = -m_kappa * r2;
423  else m_beta = 1;
424  } else {
425  m_kappa = 0;
426  if (h11 != 0 && h22 != 0) {
427  m_beta = 1 / std::sqrt(1 + h12 * h12 / h11 / h11);
428  m_alpha = std::sqrt(1 - m_beta * m_beta);
429  } else if (h11 != 0) {
430  m_beta = 1;
431  m_alpha = 0;
432  } else {
433  m_beta = 0;
434  m_alpha = 1;
435  }
436  }
437  if ((m_alpha * m_xav + m_beta * m_yav) *
438  (m_beta * m_xav - m_alpha * m_yav) < 0) neg();
439 // if (std::fabs(m_alpha)<0.01 && std::fabs(m_beta)<0.01) {
440 // std::cout << " lambda=" << lambda << " " << *this << std::endl;
441 // }
442  if (lambda >= 0) m_chisq = lambda * m_wsum_temp * m_rscale * m_rscale;
443  return lambda;
444  }
445 
446  double TRGCDCLpav::fit(double x, double y, double w)
447  {
448  if (m_nc <= 3) return -1;
449  m_chisq = -1;
450  if (m_nc < 4) {
451  calculate_average3(x, y, w);
452  double q = calculate_lpar3();
453  if (q > 0) m_chisq = q * m_wsum_temp * m_rscale * m_rscale;
454  } else {
455  calculate_average(x, y, w);
456  double q = calculate_lpar();
457  if (q > 0) m_chisq = q * m_wsum_temp * m_rscale * m_rscale;
458  }
459  return m_chisq;
460  }
461 
462  double TRGCDCLpav::fit(void)
463  {
464  if (m_nc <= 3) return -1;
465  m_chisq = -1;
466  double q;
467  if (m_nc < 4) {
468  calculate_average3();
469  q = calculate_lpar3();
470  if (q > 0) m_chisq = q * m_wsum_temp * m_rscale * m_rscale;
471  } else {
472  calculate_average();
473  q = calculate_lpar();
474  if (q > 0) m_chisq = q * m_wsum_temp * m_rscale * m_rscale;
475  }
476  return m_chisq;
477  }
478 
479  CLHEP::HepSymMatrix TRGCDCLpav::cov(int inv) const
480 #ifdef BELLE_OPTIMIZED_RETURN
481  return vret(4);
482  {
483 #else
484  {
485  CLHEP::HepSymMatrix vret(4);
486 #endif
487  vret(1, 1) = m_xxsum;
488  vret(2, 1) = m_xysum;
489  vret(2, 2) = m_yysum;
490  vret(3, 1) = m_xsum;
491  vret(3, 2) = m_ysum;
492  vret(3, 3) = m_wsum;
493  vret(4, 1) = m_xrrsum;
494  vret(4, 2) = m_yrrsum;
495  vret(4, 3) = m_xxsum + m_yysum;
496  vret(4, 4) = m_rrrrsum;
497  if (inv == 0)
498  {
499 // int i=vret.Inv();
500  int i;
501  vret.invert(i);
502  if (i != 0) {
503  std::cout << "TRGCDCLpav::cov:could not invert nc=" << m_nc << vret;
504 #ifdef HAVE_EXCEPTION
505  throw new Singular();
506 #endif
507  }
508  }
509  return vret;
510  }
511 
512  CLHEP::HepSymMatrix TRGCDCLpav::cov_c(int inv) const
513 #ifdef BELLE_OPTIMIZED_RETURN
514  return vret(3);
515  {
516 #else
517  {
518  CLHEP::HepSymMatrix vret(3);
519 #endif
520 #ifdef HAVE_EXCEPTION
521  try {
522 #endif
523  vret = cov(1).similarity(dldc());
524 #ifdef HAVE_EXCEPTION
525  } catch (TRGCDCLpav::Singular)
526  {
527  throw new Singular_c();
528  }
529 #endif
530  if (inv == 0)
531  {
532 // int i = vret.Inv();
533  int i;
534  vret.invert(i);
535  if (i != 0) {
536  std::cout << "TRGCDCLpav::cov_c:could not invert " << vret;
537 #ifdef HAVE_EXCEPTION
538  throw new Singular_c();
539 #endif
540  }
541  }
542  return vret;
543  }
544 
545  int TRGCDCLpav::extrapolate(double r, double& phi, double& dphi) const
546  {
547  double x, y;
548  if (m_chisq < 0) return -1;
549  if (xy(r, x, y) != 0) return -1;
550  phi = std::atan2(y, x);
551  if (phi < 0) phi += (2 * M_PI);
552  CLHEP::HepVector v(4);
553  v(1) = x;
554  v(2) = y;
555  v(3) = 1;
556  v(4) = r * r;
557 // CLHEP::HepSymMatrix l = cov().similarityT(v);
558 #ifdef HAVE_EXCEPTION
559  try {
560 #endif
561 // CLHEP::HepSymMatrix l = cov().similarity(v.T());
562 // // std::cout << "delta d^2=" << l(1,1);
563 // if (l(1,1)>0) {
564  double l = cov().similarity(v);
565  if (l > 0) {
566  double ls = std::sqrt(l);
567  dphi = ls / r;
568  // std::cout << " delta d=" << ls << " dphi=" << dphi;
569  }
570 #ifdef HAVE_EXCEPTION
571  } catch (TRGCDCLpav::Singular) {
572  return -1;
573  }
574 #endif
575 // std::cout << std::endl;
576  return 0;
577  }
578 
579  double TRGCDCLpav::similarity(double x, double y) const
580  {
581  if (m_nc <= 3) return -1;
582  CLHEP::HepVector v(4);
583  v(1) = x;
584  v(2) = y;
585  v(3) = 1;
586  v(4) = x * x + y * y;
587  double l;
588 #ifdef HAVE_EXCEPTION
589  try {
590 #endif
591  l = cov().similarity(v);
592 #ifdef HAVE_EXCEPTION
593  } catch (TRGCDCLpav::Singular) {
594  return -1;
595  }
596 #endif
597  return l;
598  }
599 
600 // void TRGCDCLpav::add(double xi, double yi, double w, double a, double b)
601 // {
602 //double wi = err_dis_inv(xi, yi, w, a, b);
603 // add(xi, yi, wi); calling itself with output
604 // }
605 //
606  void TRGCDCLpav::add_point(double xi, double yi,
607  double wi)
608  {
609  m_wsum += wi;
610  m_xsum += wi * xi;
611  m_ysum += wi * yi;
612  m_xxsum += wi * xi * xi;
613  m_yysum += wi * yi * yi;
614  m_xysum += wi * xi * yi;
615 //double rri = ( xi * xi + yi * yi );
616 //double wrri = wi * rri;
617  double rri = (xi * xi + yi * yi);
618  double wrri = wi * rri;
619  m_xrrsum += wrri * xi;
620  m_yrrsum += wrri * yi;
621  m_rrrrsum += wrri * rri;
622  m_nc += 1;
623  }
624 
625  void TRGCDCLpav::add_point_frac(double xi, double yi, double w, double a)
626  {
627 //double wi = w * a;
628  double wi = w * a;
629  m_wsum += wi;
630  m_xsum += wi * xi;
631  m_ysum += wi * yi;
632  m_xxsum += wi * xi * xi;
633  m_yysum += wi * yi * yi;
634  m_xysum += wi * xi * yi;
635 //double rri = ( xi * xi + yi * yi );
636 //double wrri = wi * rri;
637  double rri = (xi * xi + yi * yi);
638  double wrri = wi * rri;
639  m_xrrsum += wrri * xi;
640  m_yrrsum += wrri * yi;
641  m_rrrrsum += wrri * rri;
642  m_nc += a;
643  }
644 
645  void TRGCDCLpav::sub(double xi, double yi, double w, double a, double b)
646  {
647 //double wi = err_dis_inv(xi, yi, w, a, b);
648  double wi = err_dis_inv(xi, yi, w, a, b);
649  m_wsum -= wi;
650  m_xsum -= wi * xi;
651  m_ysum -= wi * yi;
652  m_xxsum -= wi * xi * xi;
653  m_yysum -= wi * yi * yi;
654  m_xysum -= wi * xi * yi;
655 //double rri = ( xi * xi + yi * yi );
656 //double wrri = wi * rri;
657  double rri = (xi * xi + yi * yi);
658  double wrri = wi * rri;
659  m_xrrsum -= wrri * xi;
660  m_yrrsum -= wrri * yi;
661  m_rrrrsum -= wrri * rri;
662  m_nc -= 1;
663  }
664 
665  const TRGCDCLpav& TRGCDCLpav::operator+=(const TRGCDCLpav& la1)
666  {
667  m_wsum += la1.m_wsum;
668  m_xsum += la1.m_xsum;
669  m_ysum += la1.m_ysum;
670  m_xxsum += la1.m_xxsum;
671  m_yysum += la1.m_yysum;
672  m_xysum += la1.m_xysum;
673  m_xrrsum += la1.m_xrrsum;
674  m_yrrsum += la1.m_yrrsum;
675  m_rrrrsum += la1.m_rrrrsum;
676  m_nc += la1.m_nc;
677  return *this;
678  }
679 
680  TRGCDCLpav operator+(const TRGCDCLpav& la1, const TRGCDCLpav& la2)
681 #ifdef BELLE_OPTIMIZED_RETURN
682  return la;
683  {
684 #else
685  {
686  TRGCDCLpav la;
687 #endif
688  la.m_wsum = la1.m_wsum + la2.m_wsum;
689  la.m_xsum = la1.m_xsum + la2.m_xsum;
690  la.m_ysum = la1.m_ysum + la2.m_ysum;
691  la.m_xxsum = la1.m_xxsum + la2.m_xxsum;
692  la.m_yysum = la1.m_yysum + la2.m_yysum;
693  la.m_xysum = la1.m_xysum + la2.m_xysum;
694  la.m_xrrsum = la1.m_xrrsum + la2.m_xrrsum;
695  la.m_yrrsum = la1.m_yrrsum + la2.m_yrrsum;
696  la.m_rrrrsum = la1.m_rrrrsum + la2.m_rrrrsum;
697  la.m_nc = la1.m_nc + la2.m_nc;
698  return la;
699  }
700 
701  double TRGCDCLpav::prob() const
702  {
703  if (m_nc <= 3) return 0;
704  if (m_chisq < 0) return 0;
705 //cnv float c = m_chisq;
706 //cnv int nci = (int)m_nc - 3;
707  // temporarily commented out
708  // double p = (double) prob_(&c, &nci);
709  double p = 0;
710  return p;
711  }
712 
713  double TRGCDCLpav::chi_deg() const
714  {
715  if (m_nc <= 3) return -1;
716  else return m_chisq / (m_nc - 3);
717  }
718 
719  double TRGCDCLpav::delta_chisq(double x, double y, double w) const
720  {
721  double sim = similarity(x, y);
722  if (sim < 0) return -1;
723  double d = d0(x, y);
724  double delta = sqr(d) * w / (1 + sim * w);
725  return delta;
726  }
727 
729 } // namespace Belle2
Belle2::TRGCDCLpav::m_rscale
double m_rscale
data members
Definition: Lpav.h:164
Belle2::TRGCDCLpav::m_xysum
double m_xysum
data members
Definition: Lpav.h:146
Belle2::TRGCDCLpar::move
void move(double x, double y)
private member functions
Definition: Lpar.h:231
Belle2::operator<<
std::ostream & operator<<(std::ostream &output, const IntervalOfValidity &iov)
Definition: IntervalOfValidity.cc:196
Belle2::TRGCDCLpav::m_yyavp
double m_yyavp
data members
Definition: Lpav.h:168
Belle2::TRGCDCLpar::xy
bool xy(double, double &, double &, int dir=0) const
private const member functions
Definition: Lpar.cc:162
Belle2::TRGCDCLpav::solve_lambda3
double solve_lambda3(void)
private member function solve_lambda3
Definition: Lpav.cc:324
Belle2::TRGCDCLpar::y
double y(double r) const
private const member functions
Definition: Lpar.cc:186
Belle2::TRGCDCLpav::solve_lambda
double solve_lambda(void)
private member function solve_lambda
Definition: Lpav.cc:287
Belle2::TRGCDCLpav::calculate_average
void calculate_average(void)
member functions for calculation
Definition: Lpav.cc:132
Belle2::TRGCDCLpav::m_rrrrsum
double m_rrrrsum
data members
Definition: Lpav.h:152
Belle2::TRGCDCLpav::m_xrravp
double m_xrravp
data members
Definition: Lpav.h:170
Belle2::TRGCDCLpav::fit
double fit()
member functions for fit
Definition: Lpav.cc:462
Belle2::err_dis_inv
static double err_dis_inv(double x, double y, double w, double a, double b)
distance error
Definition: Lpav.cc:35
Belle2::TRGCDCLpav::m_yav
double m_yav
data members
Definition: Lpav.h:159
Belle2::TRGCDCLpav::m_yrravp
double m_yrravp
data members
Definition: Lpav.h:172
Belle2::TRGCDCLpar::m_alpha
double m_alpha
data members
Definition: Lpar.h:175
Belle2::TRGCDCLpav::TRGCDCLpav
TRGCDCLpav()
Constructor.
Definition: Lpav.cc:54
Belle2::TRGCDCLpav::prob
double prob() const
const member function prob
Definition: Lpav.cc:701
Belle2::TRGCDCLpav::m_wsum
double m_wsum
data members
Definition: Lpav.h:136
Belle2::TRGCDCLpar
TRGCDCLpar class.
Definition: Lpar.h:36
Belle2::TRGCDCLpav::m_yrrsum
double m_yrrsum
data members
Definition: Lpav.h:150
Belle2::TRGCDCLpav::extrapolate
int extrapolate(double, double &, double &) const
const member function for extrapolation
Definition: Lpav.cc:545
Belle2::TRGCDCLpar::phi
double phi(double r, int dir=0) const
const member functions
Definition: Lpar.cc:193
Belle2::TRGCDCLpar::x
double x(double r) const
private const member functions
Definition: Lpar.cc:179
Belle2::TRGCDCLpar::neg
void neg()
member functions
Definition: Lpar.h:243
Belle2::TRGCDCLpav::m_sinrot
double m_sinrot
data members
Definition: Lpav.h:176
Belle2::TRGCDCLpar::m_gamma
double m_gamma
data members
Definition: Lpar.h:179
Belle2::TRGCDCLpav::calculate_lpar
double calculate_lpar(void)
member functions for calculation
Definition: Lpav.cc:348
Belle2::TRGCDCLpar::d0
double d0(double x, double y) const
private const member functions
Definition: Lpar.h:251
Belle2::TRGCDCLpav::m_xrrsum
double m_xrrsum
data members
Definition: Lpav.h:148
Belle2::TRGCDCLpav::m_wsum_temp
double m_wsum_temp
data members
Definition: Lpav.h:155
Belle2::TRGCDCLpav::Singular
exception class, no covarience matrix.
Definition: Lpav.h:99
Belle2::TRGCDCLpav
TRGCDCLpav class.
Definition: Lpav.h:32
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TRGCDCLpav::sub
void sub(double x, double y, double w=1, double a=0, double b=0)
private member functions
Definition: Lpav.cc:645
Belle2::TRGCDCLpav::delta_chisq
double delta_chisq(double x, double y, double w=1) const
const member function for delta chisq
Definition: Lpav.cc:719
Belle2::TRGCDCLpav::calculate_average3
void calculate_average3(void)
member functions for calculation
Definition: Lpav.cc:241
Belle2::TRGCDCLpar::d
double d(double x, double y) const
const member functions
Definition: Lpar.h:256
Belle2::TRGCDCLpav::clear
void clear()
member functions for clear
Definition: Lpav.h:253
Belle2::TRGCDCLpav::m_xxavp
double m_xxavp
data members
Definition: Lpav.h:166
Belle2::TRGCDCLpav::~TRGCDCLpav
virtual ~TRGCDCLpav()
Destructor.
Definition: Lpav.cc:86
Belle2::operator+
B2Vector3< DataType > operator+(const TVector3 &a, const B2Vector3< DataType > &b)
non-memberfunction for adding a TVector3 to a B2Vector3
Definition: B2Vector3.h:535
Belle2::TRGCDCLpav::m_xxsum
double m_xxsum
data members
Definition: Lpav.h:142
Belle2::TRGCDCLpav::cov_c
CLHEP::HepSymMatrix cov_c(int=0) const
const member function cov_c
Definition: Lpav.cc:512
Belle2::TRGCDCLpav::operator+=
const TRGCDCLpav & operator+=(const TRGCDCLpav &)
assignment operator(s)
Definition: Lpav.cc:665
Belle2::TRGCDCLpav::m_cosrot
double m_cosrot
data members
Definition: Lpav.h:178
Belle2::TRGCDCLpav::m_ysum
double m_ysum
data members
Definition: Lpav.h:140
Belle2::TRGCDCLpav::add_point
void add_point(double x, double y, double w=1)
member functions to add point
Definition: Lpav.cc:606
Belle2::TRGCDCLpav::chi_deg
double chi_deg() const
const member function chi_deg
Definition: Lpav.cc:713
Belle2::TRGCDCLpav::m_xyavp
double m_xyavp
data members
Definition: Lpav.h:161
Belle2::TRGCDCLpav::calculate_average_n
void calculate_average_n(double xxav, double yyav, double xyav, double xrrav, double yrrav, double rrrrav)
private member function calculate_average_n
Definition: Lpav.cc:150
Belle2::TRGCDCLpav::cov
CLHEP::HepSymMatrix cov(int=0) const
const member function cov
Definition: Lpav.cc:479
Belle2::TRGCDCLpav::m_xav
double m_xav
data members
Definition: Lpav.h:157
Belle2::TRGCDCLpav::m_chisq
double m_chisq
data members
Definition: Lpav.h:183
Belle2::TRGCDCLpar::m_kappa
double m_kappa
data members
Definition: Lpar.h:181
Belle2::TRGCDCLpav::calculate_lpar3
double calculate_lpar3(void)
member functions for calculation
Definition: Lpav.cc:402
Belle2::TRGCDCLpav::m_nc
double m_nc
data members
Definition: Lpav.h:181
Belle2::TRGCDCLpav::m_xsum
double m_xsum
data members
Definition: Lpav.h:138
Belle2::TRGCDCLpav::m_rrrravp
double m_rrrravp
data members
Definition: Lpav.h:174
Belle2::TRGCDCLpav::m_yysum
double m_yysum
data members
Definition: Lpav.h:144
Belle2::prob_
float prob_(float *, int *)
prob function
Belle2::TRGCDCLpar::scale
void scale(double s)
private member functions
Definition: Lpar.h:133
Belle2::TRGCDCLpar::m_beta
double m_beta
data members
Definition: Lpar.h:177
Belle2::TRGCDCLpar::rotate
void rotate(double c, double s)
private member functions
Definition: Lpar.h:223
Belle2::TRGCDCLpav::add_point_frac
void add_point_frac(double x, double y, double w, double f)
member functions to add point
Definition: Lpav.cc:625
Belle2::TRGCDCLpav::similarity
double similarity(double, double) const
const member function similarity
Definition: Lpav.cc:579