L1Field.h

Go to the documentation of this file.

#ifndef L1Field_h
#define L1Field_h 1
 
#include "CbmL1Def.h"
 
#include <iostream>
using std::cout;
using std::endl;
using std::ostream;
 
class L1FieldValue
{
 
  public:
    L1FieldValue()
        : x(0)
        , y(0)
        , z(0) {};
 
    fvec x, y, z;
 
    void Combine(L1FieldValue& B, fvec w)
    {
        x += w * (B.x - x);
        y += w * (B.y - y);
        z += w * (B.z - z);
    }
 
    friend ostream& operator<<(ostream& out, L1FieldValue& B)
    {
        return out << B.x[0] << " | " << B.y[0] << " | " << B.z[0];
    };
} _fvecalignment;
 
class L1FieldSlice
{
 
  public:
    fvec cx[21], cy[21], cz[21];   // polinom coeff.
 
    L1FieldSlice()
    {
        for (int i = 0; i < 21; i++)
            cx[i] = cy[i] = cz[i] = 0;
    }
 
    void GetFieldValue(const fvec& x, const fvec& y, L1FieldValue& B) const
    {
        fvec x2 = x * x;
        fvec y2 = y * y;
        fvec xy = x * y;
 
        fvec x3 = x2 * x;
        fvec y3 = y2 * y;
        fvec xy2 = x * y2;
        fvec x2y = x2 * y;
 
        fvec x4 = x3 * x;
        fvec y4 = y3 * y;
        fvec xy3 = x * y3;
        fvec x2y2 = x2 * y2;
        fvec x3y = x3 * y;
 
        fvec x5 = x4 * x;
        fvec y5 = y4 * y;
        fvec xy4 = x * y4;
        fvec x2y3 = x2 * y3;
        fvec x3y2 = x3 * y2;
        fvec x4y = x4 * y;
 
        B.x = cx[0] + cx[1] * x + cx[2] * y + cx[3] * x2 + cx[4] * xy + cx[5] * y2 + cx[6] * x3 + cx[7] * x2y
              + cx[8] * xy2 + cx[9] * y3 + cx[10] * x4 + cx[11] * x3y + cx[12] * x2y2 + cx[13] * xy3 + cx[14] * y4
              + cx[15] * x5 + cx[16] * x4y + cx[17] * x3y2 + cx[18] * x2y3 + cx[19] * xy4 + cx[20] * y5;
 
        B.y = cy[0] + cy[1] * x + cy[2] * y + cy[3] * x2 + cy[4] * xy + cy[5] * y2 + cy[6] * x3 + cy[7] * x2y
              + cy[8] * xy2 + cy[9] * y3 + cy[10] * x4 + cy[11] * x3y + cy[12] * x2y2 + cy[13] * xy3 + cy[14] * y4
              + cy[15] * x5 + cy[16] * x4y + cy[17] * x3y2 + cy[18] * x2y3 + cy[19] * xy4 + cy[20] * y5;
 
        B.z = cz[0] + cz[1] * x + cz[2] * y + cz[3] * x2 + cz[4] * xy + cz[5] * y2 + cz[6] * x3 + cz[7] * x2y
              + cz[8] * xy2 + cz[9] * y3 + cz[10] * x4 + cz[11] * x3y + cz[12] * x2y2 + cz[13] * xy3 + cz[14] * y4
              + cz[15] * x5 + cz[16] * x4y + cz[17] * x3y2 + cz[18] * x2y3 + cz[19] * xy4 + cz[20] * y5;
    }
} _fvecalignment;
 
class L1FieldRegion
{
 
  public:
    L1FieldRegion()
        : cx0(0)
        , cx1(0)
        , cx2(0)
        , cy0(0)
        , cy1(0)
        , cy2(0)
        , cz0(0)
        , cz1(0)
        , cz2(0)
        , z0(0)
    {}
 
    L1FieldRegion(float reg[10])
        : cx0(reg[0])
        , cx1(reg[1])
        , cx2(reg[2])
        , cy0(reg[3])
        , cy1(reg[4])
        , cy2(reg[5])
        , cz0(reg[6])
        , cz1(reg[7])
        , cz2(reg[8])
        , z0(reg[9])
    {}
 
    fvec cx0, cx1, cx2;   // Bx(z) = cx0 + cx1*(z-z0) + cx2*(z-z0)^2
    fvec cy0, cy1, cy2;   // By(z) = cy0 + cy1*(z-z0) + cy2*(z-z0)^2
    fvec cz0, cz1, cz2;   // Bz(z) = cz0 + cz1*(z-z0) + cz2*(z-z0)^2
    fvec z0;
 
    L1FieldValue Get(const fvec z)
    {
        fvec dz = (z - z0);
        fvec dz2 = dz * dz;
        L1FieldValue B;
        B.x = cx0 + cx1 * dz + cx2 * dz2;
        B.y = cy0 + cy1 * dz + cy2 * dz2;
        B.z = cz0 + cz1 * dz + cz2 * dz2;
        return B;
    }
 
    void Set(const L1FieldValue& B0,
             const fvec B0z,
             const L1FieldValue& B1,
             const fvec B1z,
             const L1FieldValue& B2,
             const fvec B2z)
    {
        z0 = B0z;
        fvec dz1 = B1z - B0z, dz2 = B2z - B0z;
        fvec det = rcp(fvec(dz1 * dz2 * (dz2 - dz1)));
        fvec w21 = -dz2 * det;
        fvec w22 = dz1 * det;
        fvec w11 = -dz2 * w21;
        fvec w12 = -dz1 * w22;
 
        fvec dB1 = B1.x - B0.x;
        fvec dB2 = B2.x - B0.x;
        cx0 = B0.x;
        cx1 = dB1 * w11 + dB2 * w12;
        cx2 = dB1 * w21 + dB2 * w22;
 
        dB1 = B1.y - B0.y;
        dB2 = B2.y - B0.y;
        cy0 = B0.y;
        cy1 = dB1 * w11 + dB2 * w12;
        cy2 = dB1 * w21 + dB2 * w22;
 
        dB1 = B1.z - B0.z;
        dB2 = B2.z - B0.z;
        cz0 = B0.z;
        cz1 = dB1 * w11 + dB2 * w12;
        cz2 = dB1 * w21 + dB2 * w22;
    }
 
    void Set(const L1FieldValue& B0, const fvec B0z, const L1FieldValue& B1, const fvec B1z)
    {
        z0 = B0z[0];
        fvec dzi = rcp(fvec(B1z - B0z));
        cx0 = B0.x;
        cy0 = B0.y;
        cz0 = B0.z;
        cx1 = (B1.x - B0.x) * dzi;
        cy1 = (B1.y - B0.y) * dzi;
        cz1 = (B1.z - B0.z) * dzi;
        cx2 = cy2 = cz2 = 0;
    }
 
    void Shift(fvec z)
    {
        fvec dz = z - z0;
        fvec cx2dz = cx2 * dz;
        fvec cy2dz = cy2 * dz;
        fvec cz2dz = cz2 * dz;
        z0 = z[0];
        cx0 += (cx1 + cx2dz) * dz;
        cy0 += (cy1 + cy2dz) * dz;
        cz0 += (cz1 + cz2dz) * dz;
        cx1 += cx2dz + cx2dz;
        cy1 += cy2dz + cy2dz;
        cz1 += cz2dz + cz2dz;
    }
 
    void SetOneEntry(const int i0, const L1FieldRegion& f1, const int i1)
    {
        cx0[i0] = f1.cx0[i1];
        cx1[i0] = f1.cx1[i1];
        cx2[i0] = f1.cx2[i1];
        cy0[i0] = f1.cy0[i1];
        cy1[i0] = f1.cy1[i1];
        cy2[i0] = f1.cy2[i1];
        cz0[i0] = f1.cz0[i1];
        cz1[i0] = f1.cz1[i1];
        cz2[i0] = f1.cz2[i1];
        z0[i0] = f1.z0[i1];
    }
 
    void SetOneEntry(const L1FieldRegion& f1, const int i1)
    {
        cx0 = f1.cx0[i1];
        cx1 = f1.cx1[i1];
        cx2 = f1.cx2[i1];
        cy0 = f1.cy0[i1];
        cy1 = f1.cy1[i1];
        cy2 = f1.cy2[i1];
        cz0 = f1.cz0[i1];
        cz1 = f1.cz1[i1];
        cz2 = f1.cz2[i1];
        z0 = f1.z0[i1];
    }
 
    void GetOneEntry(float reg[10], const int iVec)
    {
        reg[0] = cx0[iVec];
        reg[1] = cx1[iVec];
        reg[2] = cx2[iVec];
        reg[3] = cy0[iVec];
        reg[4] = cy1[iVec];
        reg[5] = cy2[iVec];
        reg[6] = cz0[iVec];
        reg[7] = cz1[iVec];
        reg[8] = cz2[iVec];
        reg[9] = z0[iVec];
    }
} _fvecalignment;
 
#endif