Changeset 4982 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
09/13/04 15:34:38 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mbase/MMath.cc

    r4966 r4982  
    126126// If the determinant==0 an empty TVector3 is returned.
    127127//
     128#include <iostream>
    128129TVector3 MMath::GetParab(const TVector3 &x, const TVector3 &y)
    129130{
    130     const TVector3 sq(x(0)*x(0), x(1)*x(1), x(2)*x(2));
     131    Double_t x1 = x(0);
     132    Double_t x2 = x(1);
     133    Double_t x3 = x(2);
    131134
    132     const TVector3 ai2 = sq.Cross(sq);
     135    Double_t y1 = y(0);
     136    Double_t y2 = y(1);
     137    Double_t y3 = y(2);
    133138
    134     const Double_t det = x.Dot(ai2);
     139    const double det =
     140        + x2*x3*x3 + x1*x2*x2 + x3*x1*x1
     141        - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;
     142
     143
    135144    if (det==0)
    136145        return TVector3();
    137146
    138     const TVector3 ai1 = x.Cross(sq);
    139     const TVector3 ai3 = x.Cross(x);
     147    const double det1 = 1.0/det;
    140148
    141     TVector3 res(y.Dot(ai1), y.Dot(ai2), y.Dot(ai3));
    142     res *= 1./det;
     149    const double ai11 = x2*x3*x3 - x3*x2*x2;
     150    const double ai12 = x3*x1*x1 - x1*x3*x3;
     151    const double ai13 = x1*x2*x2 - x2*x1*x1;
    143152
    144     return res;
     153    const double ai21 = x2*x2 - x3*x3;
     154    const double ai22 = x3*x3 - x1*x1;
     155    const double ai23 = x1*x1 - x2*x2;
     156
     157    const double ai31 = x3 - x2;
     158    const double ai32 = x1 - x3;
     159    const double ai33 = x2 - x1;
     160
     161    return TVector3((ai11*y1 + ai12*y2 + ai13*y3) * det1,
     162                    (ai21*y1 + ai22*y2 + ai23*y3) * det1,
     163                    (ai31*y1 + ai32*y2 + ai33*y3) * det1);
    145164}
    146165
     
    158177
    159178    const TVector3 vx0(l0, l1, l2);
    160     return pow(10, InterpolParabLin(vx0, vy, TMath::Log10(x)));
     179    return InterpolParabLin(vx0, vy, TMath::Log10(x));
    161180}
    162181
     
    168187
    169188    const TVector3 vx0(l0, l1, l2);
    170     return TMath::ACos(InterpolParabLin(vx0, vy, TMath::Cos(x)));
     189    return InterpolParabLin(vx0, vy, TMath::Cos(x));
    171190}
Note: See TracChangeset for help on using the changeset viewer.