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