FEDRA emulsion software from the OPERA Collaboration
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | List of all members
EdbMath Class Reference

general matematical algorithms More...

#include <EdbMath.h>

Public Member Functions

 EdbMath ()
 
virtual ~EdbMath ()
 

Static Public Member Functions

static double Angle3 (float tx1, float ty1, float tx2, float ty2)
 
static int ArrStat (int n, float *x, float par[4])
 
static double DistancePointLine3 (float Point[3], float LineStart[3], float LineEnd[3], bool inside)
 
static int LFIT3 (float *X, float *Y, float *Z, float *W, int L, float &X0, float &Y0, float &Z0, float &TX, float &TY, float &EX, float &EY)
 
static void LFITW (float *X, float *Y, float *W, int L, int KEY, float &A, float &B, float &E)
 
static bool LineLineIntersect (float p1[3], float p2[3], float p3[3], float p4[3], float pa[3], float pb[3], double &mua, double &mub)
 
static bool LinFitCOV (int n, float *x, float *y, double *c, float *p, float *d, float *chi2)
 
static bool LinFitDiag (int n, float *x, float *y, float *e, float p[2], float d[2][2], float *chi2)
 
static double Magnitude3 (float Point1[3], float Point2[3])
 

Detailed Description

general matematical algorithms

//////////////////////////////////////////////////////////////////////// // EdbMath // // Collection of matematical algorithns // // ////////////////////////////////////////////////////////////////////////

Constructor & Destructor Documentation

◆ EdbMath()

EdbMath::EdbMath ( )
inline
21{}

◆ ~EdbMath()

virtual EdbMath::~EdbMath ( )
inlinevirtual
22{}

Member Function Documentation

◆ Angle3()

double EdbMath::Angle3 ( float  tx1,
float  ty1,
float  tx2,
float  ty2 
)
static

calculate space angle between two vectors defined by the direction tangents

43{
45 double a2 = tx1*tx1+ty1*ty1+1;
46 double b2 = tx2*tx2+ty2*ty2+1;
47 double c2 = tx1*tx2+ty1*ty2+1;
48 return ACos( c2/Sqrt(a2*b2) );
49}
c2
TCanvas * c2
Definition: energy.C:26

◆ ArrStat()

int EdbMath::ArrStat ( int  n,
float *  x,
float  par[4] 
)
static

calculate array statistics:
output: par (mean, rms, xmin, xmax);

24{
27
28 if(n<1) return 0;
29 par[0]=par[2]=par[3]=x[0];
30 double xs=0; for(int i=0; i<n; i++) xs+=x[i];
31 double x0 = xs/n;
32 double dxs=0; for(int i=0; i<n; i++) dxs += ((x[i]-x0)*(x[i]-x0));
33 double rms = Sqrt(dxs/n);
34 for(int i=0; i<n; i++) par[2] = par[2]<x[i]?par[2]:x[i];
35 for(int i=0; i<n; i++) par[3] = par[3]>x[i]?par[3]:x[i];
36 par[0]=x0;
37 par[1]=rms;
38 return n;
39}

◆ DistancePointLine3()

double EdbMath::DistancePointLine3 ( float  Point[3],
float  LineStart[3],
float  LineEnd[3],
bool  inside 
)
static

Disclamer: This is a modified code from site of Paul Bourke

62{
64
65 double LineMag;
66 double U;
67 float Intersection[3];
68
69 LineMag = Magnitude3( LineEnd, LineStart );
70
71 U = ( ( ( Point[0] - LineStart[0] ) * ( LineEnd[0] - LineStart[0] ) ) +
72 ( ( Point[1] - LineStart[1] ) * ( LineEnd[1] - LineStart[1] ) ) +
73 ( ( Point[2] - LineStart[2] ) * ( LineEnd[2] - LineStart[2] ) ) ) /
74 ( LineMag * LineMag );
75
76 if( U < 0.0f || U > 1.0f )
77 inside=false; // closest point does not fall within the line segment
78 else inside=true;
79
80 Intersection[0] = LineStart[0] + U * ( LineEnd[0] - LineStart[0] );
81 Intersection[1] = LineStart[1] + U * ( LineEnd[1] - LineStart[1] );
82 Intersection[2] = LineStart[2] + U * ( LineEnd[2] - LineStart[2] );
83
84 return Magnitude3( Point, Intersection );
85}
EdbMath::Magnitude3
static double Magnitude3(float Point1[3], float Point2[3])
Definition: EdbMath.cxx:52

◆ LFIT3()

int EdbMath::LFIT3 ( float *  X,
float *  Y,
float *  Z,
float *  W,
int  L,
float &  X0,
float &  Y0,
float &  Z0,
float &  TX,
float &  TY,
float &  EX,
float &  EY 
)
static

Linar fit in 3-d case (microtrack-like)
Input: X,Y,Z - coords, W -weight - arrays of the lengh >=L
Note that X,Y,Z modified by function
Output: X0,Y0,Z0 - center of gravity of segment
TX,TY : tangents in respect to Z axis

315{
321
322 if(L<2) return 0;
323
324 // first find C.O.G.:
325 double wx=0, wy=0, wz=0, w=0;
326 for(int i=0; i<L; i++) {
327 wx += X[i]*W[i];
328 wy += Y[i]*W[i];
329 wz += Z[i]*W[i];
330 w += W[i];
331 }
332 X0 = wx/w;
333 Y0 = wy/w;
334 Z0 = wz/w;
335 for(int i=0; i<L; i++) {
336 X[i] -= X0;
337 Y[i] -= Y0;
338 Z[i] -= Z0;
339 }
340
341 if(L==2) {
342 TX = (X[1]-X[0])/(Z[1]-Z[0]);
343 TY = (Y[1]-Y[0])/(Z[1]-Z[0]);
344 EX=EY=0;
345 return L;
346 }
347
348 float x0=0,y0=0;
349 EdbMath::LFITW( Z,X,W, L, 1, TX,x0,EX );
350 EdbMath::LFITW( Z,Y,W, L, 1, TY,y0,EY );
351
352 X0 += x0;
353 Y0 += y0;
354
355 return L;
356}
X0
brick X0
Definition: RecDispMC.C:112
EdbMath::LFITW
static void LFITW(float *X, float *Y, float *W, int L, int KEY, float &A, float &B, float &E)
Definition: EdbMath.cxx:262
Z0
float Z0
Definition: hwinit.C:67
Header::Center::X
Double_t X
Definition: tlg2couples.C:76
Header::Center::Y
Double_t Y
Definition: tlg2couples.C:76
Header::Transform::TY
Double_t TY
Definition: tlg2couples.C:78
Header::Transform::TX
Double_t TX
Definition: tlg2couples.C:78
MIPEmulsionTrack::Z
Double_t Z
Definition: tlg2couples.C:104
W
Int_t W
Definition: testBGReduction_By_ANN.C:15
w
void w(int rid=2, int nviews=2)
Definition: test.C:27

◆ LFITW()

void EdbMath::LFITW ( float *  X,
float *  Y,
float *  W,
int  L,
int  KEY,
float &  A,
float &  B,
float &  E 
)
static

E250:
TO PERFORM A WEIGHTED STRAIGHT LINE FIT
FOR FORMULAE USED SEE MENZEL, FORMULAS OF PHYSICS P.116
FIT IS OF Y=AX+B , WITH S**2 ESTIMATOR E. WEIGHTS ARE IN W.
IF KEY=0, POINTS WITH Y=0 ARE IGNORED
L IS NO. OF POINTS

263{
273
274 A=0; B=0; E=0;
275
276 // CALCULTE SUMS
277 if(L<2) return;
278
279 double WW=0.;
280 double WWF=0.;
281 double WWFI=0.;
282 double W2=0.;
283 double W2X=0.;
284 double W2Y=0.;
285 double W2XY=0.;
286 double W2X2=0.;
287 double W2Y2=0.;
288 int ICNT=0;
289
290 for( int j=0; j<L; j++ ){
291 if( Y[j]==0. && KEY==0 ) continue;
292 WW=W[j]*W[j];
293 W2=WW+W2;
294 WWF=WW*X[j];
295 W2X=WWF+W2X;
296 W2X2=WWF*X[j]+W2X2;
297 W2XY=WWF*Y[j]+W2XY;
298 WWFI=WW*Y[j];
299 W2Y=WWFI+W2Y;
300 W2Y2=WWFI*Y[j]+W2Y2;
301 ICNT++;
302 }
303
304 // FIT PARAMETERS
305 A=(W2XY-W2X*W2Y/W2)/(W2X2-W2X*W2X/W2);
306 B=(W2Y-A*W2X)/W2;
307 if(ICNT<3) return;
308 E = (W2Y2-W2Y*W2Y/W2 -
309 (W2XY-W2X*W2Y/W2)*(W2XY-W2X*W2Y/W2)/(W2X2-W2X*W2X/W2))/(ICNT-2);
310}

◆ LineLineIntersect()

bool EdbMath::LineLineIntersect ( float  p1[3],
float  p2[3],
float  p3[3],
float  p4[3],
float  pa[3],
float  pb[3],
double &  mua,
double &  mub 
)
static

Disclamer: This is a modified code from site of Paul Bourke 

Calculate the line segment PaPb that is the shortest route between \n
two lines P1P2 and P3P4. Calculate also the values of mua and mub where \n
   Pa = P1 + mua (P2 - P1) \n
   Pb = P3 + mub (P4 - P3) \n
Return FALSE if no solution exists (the lines are parallel)
91{
99
100 float p13[3],p43[3],p21[3];
101 double d1343,d4321,d1321,d4343,d2121;
102 double numer,denom;
103 const float EPS=1.E-6;
104
105 p13[0] = p1[0] - p3[0];
106 p13[1] = p1[1] - p3[1];
107 p13[2] = p1[2] - p3[2];
108 p43[0] = p4[0] - p3[0];
109 p43[1] = p4[1] - p3[1];
110 p43[2] = p4[2] - p3[2];
111 if (TMath::Abs(p43[0]) < EPS &&
112 TMath::Abs(p43[1]) < EPS &&
113 TMath::Abs(p43[2]) < EPS) return false;
114 p21[0] = p2[0] - p1[0];
115 p21[1] = p2[1] - p1[1];
116 p21[2] = p2[2] - p1[2];
117 if (TMath::Abs(p21[0]) < EPS &&
118 TMath::Abs(p21[1]) < EPS &&
119 TMath::Abs(p21[2]) < EPS) return false;
120
121 d1343 = p13[0] * p43[0] + p13[1] * p43[1] + p13[2] * p43[2];// e
122 d4321 = p43[0] * p21[0] + p43[1] * p21[1] + p43[2] * p21[2];// b
123 d1321 = p13[0] * p21[0] + p13[1] * p21[1] + p13[2] * p21[2];// d
124 d4343 = p43[0] * p43[0] + p43[1] * p43[1] + p43[2] * p43[2];// c
125 d2121 = p21[0] * p21[0] + p21[1] * p21[1] + p21[2] * p21[2];// a
126
127 // a c b b
128 denom = d2121 * d4343 - d4321 * d4321;
129 if (TMath::Abs(denom) < EPS) return false;
130
131 numer = d1343 * d4321 - d1321 * d4343;
132
133 mua = numer / denom;
134 mub = (d1343 + d4321 * (mua)) / d4343;
135
136 pa[0] = p1[0] + mua * p21[0];
137 pa[1] = p1[1] + mua * p21[1];
138 pa[2] = p1[2] + mua * p21[2];
139 pb[0] = p3[0] + mub * p43[0];
140 pb[1] = p3[1] + mub * p43[1];
141 pb[2] = p3[2] + mub * p43[2];
142
143 return true;
144}

◆ LinFitCOV()

bool EdbMath::LinFitCOV ( int  n,
float *  x,
float *  y,
double *  c,
float *  p,
float *  d,
float *  chi2 
)
static

linear fit by YP with using full covariance matrix of measurements :

y = p[0] + p[1]*x , d[2][2] is error matrix for p[2]

Input : n - number of points

x,y - data arrays of the length n

c - full covariance matrix c[n][n]

194{
200
201 if (n < 2) return false;
202 int i = 0, j = 0;
203 double dy = 0;
204
205 TMatrixD D(n,n); // error matrix of measurements
206 TMatrixD A(n,2); // construction matrix
207 TMatrixD AT(2,n); // A transposed
208 TVectorD Y(n); // measurements vector
209 TVectorD Y1(n); // temporary vector
210 TVectorD DY(n); // residuals
211 TMatrixD DI(n,n); // Inverse of D-matrix
212 TMatrixD DP(2,2); // error matrix of parameters
213 TMatrixD B(2,2); // inverse error matrix of parameters
214 TVectorD P(2); // parameters vector
215
216 for (i=0; i<n; i++)
217 {
218 A(i,0) = (double)1.;
219 A(i,1) = (double)x[i];
220 AT(0,i) = A(i,0);
221 AT(1,i) = A(i,1);
222 Y(i) = (double)y[i];
223 for (j=0; j<n; j++)
224 {
225 D(i,j) = (double)(*(c+i*n+j));
226 }
227 }
228
229 Y1 = Y;
230 DI = D.Invert();
231 B = AT*(DI*A);
232 DP = B.Invert();
233 Y1 *= DI;
234 AT.ResizeTo(n,n);
235 Y1 *= AT;
236 Y1.ResizeTo(2);
237 P = Y1;
238 P *= DP;
239
240 for (i=0; i<2; i++)
241 {
242 *(p+i) = (float)P(i);
243 for (j=0; j<2; j++)
244 {
245 *(d+i*2+j) = (float)DP(i,j);
246 }
247 }
248
249 for (i=0; i<n; i++)
250 {
251 dy = (double)(y[i] - *p - (*(p+1))*x[i]);
252 DY(i) = dy;
253 }
254 TVectorD DY1 = DY;
255 DY *= DI;
256 *chi2 = (float)(DY1*DY);
257
258 return true;
259}
d
void d()
Definition: RecDispEX.C:381
p
p
Definition: testBGReduction_AllMethods.C:8
chi2
Float_t chi2
Definition: testBGReduction_By_ANN.C:14

◆ LinFitDiag()

bool EdbMath::LinFitDiag ( int  n,
float *  x,
float *  y,
float *  e,
float  p[2],
float  d[2][2],
float *  chi2 
)
static

linear fit by YP: y = p[0] + p[1]*x , d[2][2] is error matrix for p[2]
Input : n - number of points
x,y,e - data arrays of the length n

148{
152
153 if (n < 2) return false;
154 int i = 0;
155 double w = 0., ed = 0.;
156 double sw = 0., swx = 0., swx2 = 0., swy = 0., swxy = 0.;
157 float dy = 0.;
158
159 for (i=0; i<n; i++)
160 {
161 ed = (double)e[i];
162 if (ed > 0.)
163 {
164 w = 1./(ed*ed);
165 sw += w;
166 swx += w*x[i];
167 swx2 += w*x[i]*x[i];
168 swy += w*y[i];
169 swxy += w*x[i]*y[i];
170 }
171 }
172 double det = sw*swx2 - swx*swx;
173 if (det == 0.) return false;
174 d[0][0] = (float)( swx2/det);
175 d[0][1] = (float)(-swx/det);
176 d[1][0] = d[0][1];
177 d[1][1] = (float)(sw/det);
178 p[0] = (float)((swy*swx2 - swxy*swx)/det);
179 p[1] = (float)((swxy*sw - swy*swx)/det);
180 *chi2 = 0.;
181 for (i=0; i<n; i++)
182 {
183 if ( e[i] > 0.)
184 {
185 dy = y[i] - p[0] - p[1]*x[i];
186 *chi2 += dy*dy/(e[i]*e[i]);
187 }
188 }
189 return true;
190}

◆ Magnitude3()

double EdbMath::Magnitude3 ( float  Point1[3],
float  Point2[3] 
)
static
53{
54 double dx = Point2[0] - Point1[0];
55 double dy = Point2[1] - Point1[1];
56 double dz = Point2[2] - Point1[2];
57 return TMath::Sqrt( dx*dx+dy*dy+dz*dz );
58}
dz
brick dz
Definition: RecDispMC.C:107
The documentation for this class was generated from the following files: