diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index a9e5be41c11c0552a33cb44231c58da9c2ab4dd8..6b07858575b42d36875dc8402d4c656c2b27307f 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -48,9 +48,6 @@ set(SRCS
L1Algo/L1Station.cxx
L1Algo/L1TrackParFit.cxx
L1Algo/L1MCEvent.cxx
- L1Algo/L1Fit.cxx
- L1Algo/L1FitMaterial.cxx
- L1Algo/L1Extrapolation.cxx
CbmL1MCTrack.cxx
L1Algo/L1Material.cxx
L1Algo/L1UMeasurementInfo.cxx
diff --git a/reco/L1/L1Algo/L1Extrapolation.cxx b/reco/L1/L1Algo/L1Extrapolation.cxx
deleted file mode 100644
index 99ea46080f2fe13f6ef62abc1cf777cfdba8abb2..0000000000000000000000000000000000000000
--- a/reco/L1/L1Algo/L1Extrapolation.cxx
+++ /dev/null
@@ -1,916 +0,0 @@
-/* Copyright (C) 2007-2017 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Maksym Zyzak, Sergey Gorbunov [committer], Igor Kulakov */
-
-#include "L1Extrapolation.h"
-
-#include "L1Def.h"
-#include "L1Field.h"
-#include "L1TrackPar.h"
-
-
-//#define cnst static const fvec
-#define cnst const fvec
-
-void L1ExtrapolateAnalytic(L1TrackPar& T, fvec z_out, fvec qp0, const L1FieldRegion& F, const fvec* w);
-void L1ExtrapolateRungeKutta(L1TrackPar& T, fvec z_out, fvec qp0, const L1FieldRegion& F, const fvec* w);
-
-void L1Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
- (L1TrackPar& T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
- fvec z_out, // extrapolate to this z position
- fvec qp0, // use Q/p linearisation at this value
- const L1FieldRegion& F, const fvec* w)
-{
- L1ExtrapolateRungeKutta(T, z_out, qp0, F, w);
- //L1ExtrapolateAnalytic(T, z_out, qp0, F, w);
-}
-
-
-void L1ExtrapolateAnalytic(L1TrackPar& T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
- fvec z_out, // extrapolate to this z position
- fvec qp0, // use Q/p linearisation at this value
- const L1FieldRegion& F, const fvec* w)
-{
- //cout<<"Extrapolation..."<<endl;
- //
- // Part of the analytic extrapolation formula with error (c_light*B*dz)^4/4!
- //
-
- cnst c_light(0.000299792458);
-
- cnst c1(1.), c2(2.), c3(3.), c4(4.), c6(6.), c9(9.), c15(15.), c18(18.), c45(45.), c2i(1. / 2.), c3i(1. / 3.),
- c6i(1. / 6.), c12i(1. / 12.);
-
- const fvec qp = T.qp;
- fvec dz = (z_out - T.z);
-
- if (w) { dz.setZero(*w <= fvec::Zero()); }
-
- const fvec dz2 = dz * dz;
- const fvec dz3 = dz2 * dz;
-
- // construct coefficients
-
- const fvec x = T.tx;
- const fvec y = T.ty;
- const fvec xx = x * x;
- const fvec xy = x * y;
- const fvec yy = y * y;
- const fvec y2 = y * c2;
- const fvec x2 = x * c2;
- const fvec x4 = x * c4;
- const fvec xx31 = xx * c3 + c1;
- const fvec xx159 = xx * c15 + c9;
-
- const fvec Ay = -xx - c1;
- const fvec Ayy = x * (xx * c3 + c3);
- const fvec Ayz = -c2 * xy;
- const fvec Ayyy = -(c15 * xx * xx + c18 * xx + c3);
-
- const fvec Ayy_x = c3 * xx31;
- const fvec Ayyy_x = -x4 * xx159;
-
- const fvec Bx = yy + c1;
- const fvec Byy = y * xx31;
- const fvec Byz = c2 * xx + c1;
- const fvec Byyy = -xy * xx159;
-
- const fvec Byy_x = c6 * xy;
- const fvec Byyy_x = -y * (c45 * xx + c9);
- const fvec Byyy_y = -x * xx159;
-
- // end of coefficients calculation
-
- const fvec t2 = c1 + xx + yy;
- const fvec t = sqrt(t2);
- const fvec h = qp0 * c_light;
- const fvec ht = h * t;
-
- // get field integrals
- const fvec ddz = T.z - F.z0;
- const fvec Fx0 = F.cx0 + F.cx1 * ddz + F.cx2 * ddz * ddz;
- const fvec Fx1 = (F.cx1 + c2 * F.cx2 * ddz) * dz;
- const fvec Fx2 = F.cx2 * dz2;
- const fvec Fy0 = F.cy0 + F.cy1 * ddz + F.cy2 * ddz * ddz;
- const fvec Fy1 = (F.cy1 + c2 * F.cy2 * ddz) * dz;
- const fvec Fy2 = F.cy2 * dz2;
- const fvec Fz0 = F.cz0 + F.cz1 * ddz + F.cz2 * ddz * ddz;
- const fvec Fz1 = (F.cz1 + c2 * F.cz2 * ddz) * dz;
- const fvec Fz2 = F.cz2 * dz2;
-
- //
-
- const fvec sx = (Fx0 + Fx1 * c2i + Fx2 * c3i);
- const fvec sy = (Fy0 + Fy1 * c2i + Fy2 * c3i);
- const fvec sz = (Fz0 + Fz1 * c2i + Fz2 * c3i);
-
- const fvec Sx = (Fx0 * c2i + Fx1 * c6i + Fx2 * c12i);
- const fvec Sy = (Fy0 * c2i + Fy1 * c6i + Fy2 * c12i);
- const fvec Sz = (Fz0 * c2i + Fz1 * c6i + Fz2 * c12i);
-
- fvec syz;
- {
- constexpr double d = 360.;
- cnst c00(30. * 6. / d), c01(30. * 2. / d), c02(30. / d), c10(3. * 40. / d), c11(3. * 15. / d), c12(3. * 8. / d),
- c20(2. * 45. / d), c21(2. * 2. * 9. / d), c22(2. * 2. * 5. / d);
- syz = Fy0 * (c00 * Fz0 + c01 * Fz1 + c02 * Fz2) + Fy1 * (c10 * Fz0 + c11 * Fz1 + c12 * Fz2)
- + Fy2 * (c20 * Fz0 + c21 * Fz1 + c22 * Fz2);
- }
-
- fvec Syz;
- {
- constexpr double d = 2520.;
- cnst c00(21. * 20. / d), c01(21. * 5. / d), c02(21. * 2. / d), c10(7. * 30. / d), c11(7. * 9. / d),
- c12(7. * 4. / d), c20(2. * 63. / d), c21(2. * 21. / d), c22(2. * 10. / d);
- Syz = Fy0 * (c00 * Fz0 + c01 * Fz1 + c02 * Fz2) + Fy1 * (c10 * Fz0 + c11 * Fz1 + c12 * Fz2)
- + Fy2 * (c20 * Fz0 + c21 * Fz1 + c22 * Fz2);
- }
-
- const fvec syy = sy * sy * c2i;
- const fvec syyy = syy * sy * c3i;
-
- fvec Syy;
- {
- constexpr double d = 2520.;
- cnst c00(420. / d), c01(21. * 15. / d), c02(21. * 8. / d), c03(63. / d), c04(70. / d), c05(20. / d);
- Syy = Fy0 * (c00 * Fy0 + c01 * Fy1 + c02 * Fy2) + Fy1 * (c03 * Fy1 + c04 * Fy2) + c05 * Fy2 * Fy2;
- }
-
- fvec Syyy;
- {
- constexpr double d = 181440.;
- cnst c000(7560 / d), c001(9 * 1008 / d), c002(5 * 1008 / d), c011(21 * 180 / d), c012(24 * 180 / d),
- c022(7 * 180 / d), c111(540 / d), c112(945 / d), c122(560 / d), c222(112 / d);
- const fvec Fy22 = Fy2 * Fy2;
- Syyy = Fy0 * (Fy0 * (c000 * Fy0 + c001 * Fy1 + c002 * Fy2) + Fy1 * (c011 * Fy1 + c012 * Fy2) + c022 * Fy22)
- + Fy1 * (Fy1 * (c111 * Fy1 + c112 * Fy2) + c122 * Fy22) + c222 * Fy22 * Fy2;
- }
-
-
- const fvec sA1 = sx * xy + sy * Ay + sz * y;
- const fvec sA1_x = sx * y - sy * x2;
- const fvec sA1_y = sx * x + sz;
-
- const fvec sB1 = sx * Bx - sy * xy - sz * x;
- const fvec sB1_x = -sy * y - sz;
- const fvec sB1_y = sx * y2 - sy * x;
-
- const fvec SA1 = Sx * xy + Sy * Ay + Sz * y;
- const fvec SA1_x = Sx * y - Sy * x2;
- const fvec SA1_y = Sx * x + Sz;
-
- const fvec SB1 = Sx * Bx - Sy * xy - Sz * x;
- const fvec SB1_x = -Sy * y - Sz;
- const fvec SB1_y = Sx * y2 - Sy * x;
-
-
- const fvec sA2 = syy * Ayy + syz * Ayz;
- const fvec sA2_x = syy * Ayy_x - syz * y2;
- const fvec sA2_y = -syz * x2;
- const fvec sB2 = syy * Byy + syz * Byz;
- const fvec sB2_x = syy * Byy_x + syz * x4;
- const fvec sB2_y = syy * xx31;
-
- const fvec SA2 = Syy * Ayy + Syz * Ayz;
- const fvec SA2_x = Syy * Ayy_x - Syz * y2;
- const fvec SA2_y = -Syz * x2;
- const fvec SB2 = Syy * Byy + Syz * Byz;
- const fvec SB2_x = Syy * Byy_x + Syz * x4;
- const fvec SB2_y = Syy * xx31;
-
- const fvec sA3 = syyy * Ayyy;
- const fvec sA3_x = syyy * Ayyy_x;
- const fvec sB3 = syyy * Byyy;
- const fvec sB3_x = syyy * Byyy_x;
- const fvec sB3_y = syyy * Byyy_y;
-
-
- const fvec SA3 = Syyy * Ayyy;
- const fvec SA3_x = Syyy * Ayyy_x;
- const fvec SB3 = Syyy * Byyy;
- const fvec SB3_x = Syyy * Byyy_x;
- const fvec SB3_y = Syyy * Byyy_y;
-
- const fvec ht1 = ht * dz;
- const fvec ht2 = ht * ht * dz2;
- const fvec ht3 = ht * ht * ht * dz3;
- const fvec ht1sA1 = ht1 * sA1;
- const fvec ht1sB1 = ht1 * sB1;
- const fvec ht1SA1 = ht1 * SA1;
- const fvec ht1SB1 = ht1 * SB1;
- const fvec ht2sA2 = ht2 * sA2;
- const fvec ht2SA2 = ht2 * SA2;
- const fvec ht2sB2 = ht2 * sB2;
- const fvec ht2SB2 = ht2 * SB2;
- const fvec ht3sA3 = ht3 * sA3;
- const fvec ht3sB3 = ht3 * sB3;
- const fvec ht3SA3 = ht3 * SA3;
- const fvec ht3SB3 = ht3 * SB3;
-
- T.x += (x + ht1SA1 + ht2SA2 + ht3SA3) * dz;
- T.y += (y + ht1SB1 + ht2SB2 + ht3SB3) * dz;
- T.tx += ht1sA1 + ht2sA2 + ht3sA3;
- T.ty += ht1sB1 + ht2sB2 + ht3sB3;
- T.z += dz;
-
- const fvec ctdz = c_light * t * dz;
- const fvec ctdz2 = c_light * t * dz2;
-
- const fvec dqp = qp - qp0;
- const fvec t2i = c1 / t2;
- const fvec xt2i = x * t2i;
- const fvec yt2i = y * t2i;
- const fvec tmp0 = ht1SA1 + c2 * ht2SA2 + c3 * ht3SA3;
- const fvec tmp1 = ht1SB1 + c2 * ht2SB2 + c3 * ht3SB3;
- const fvec tmp2 = ht1sA1 + c2 * ht2sA2 + c3 * ht3sA3;
- const fvec tmp3 = ht1sB1 + c2 * ht2sB2 + c3 * ht3sB3;
-
- const fvec j02 = dz * (c1 + xt2i * tmp0 + ht1 * SA1_x + ht2 * SA2_x + ht3 * SA3_x);
- const fvec j12 = dz * (xt2i * tmp1 + ht1 * SB1_x + ht2 * SB2_x + ht3 * SB3_x);
- const fvec j22 = c1 + xt2i * tmp2 + ht1 * sA1_x + ht2 * sA2_x + ht3 * sA3_x;
- const fvec j32 = xt2i * tmp3 + ht1 * sB1_x + ht2 * sB2_x + ht3 * sB3_x;
-
- const fvec j03 = dz * (yt2i * tmp0 + ht1 * SA1_y + ht2 * SA2_y);
- const fvec j13 = dz * (c1 + yt2i * tmp1 + ht1 * SB1_y + ht2 * SB2_y + ht3 * SB3_y);
- const fvec j23 = yt2i * tmp2 + ht1 * sA1_y + ht2 * sA2_y;
- const fvec j33 = c1 + yt2i * tmp3 + ht1 * sB1_y + ht2 * sB2_y + ht3 * sB3_y;
-
- const fvec j04 = ctdz2 * (SA1 + c2 * ht1 * SA2 + c3 * ht2 * SA3);
- const fvec j14 = ctdz2 * (SB1 + c2 * ht1 * SB2 + c3 * ht2 * SB3);
- const fvec j24 = ctdz * (sA1 + c2 * ht1 * sA2 + c3 * ht2 * sA3);
- const fvec j34 = ctdz * (sB1 + c2 * ht1 * sB2 + c3 * ht2 * sB3);
-
-
- // extrapolate inverse momentum
-
- T.x += j04 * dqp;
- T.y += j14 * dqp;
- T.tx += j24 * dqp;
- T.ty += j34 * dqp;
-
- // covariance matrix transport
-
- const fvec c42 = T.C42, c43 = T.C43;
-
- const fvec cj00 = T.C00 + T.C20 * j02 + T.C30 * j03 + T.C40 * j04;
- // const fvec cj10 = T.C10 + T.C21*j02 + T.C31*j03 + T.C41*j04;
- const fvec cj20 = T.C20 + T.C22 * j02 + T.C32 * j03 + c42 * j04;
- const fvec cj30 = T.C30 + T.C32 * j02 + T.C33 * j03 + c43 * j04;
-
- const fvec cj01 = T.C10 + T.C20 * j12 + T.C30 * j13 + T.C40 * j14;
- const fvec cj11 = T.C11 + T.C21 * j12 + T.C31 * j13 + T.C41 * j14;
- const fvec cj21 = T.C21 + T.C22 * j12 + T.C32 * j13 + c42 * j14;
- const fvec cj31 = T.C31 + T.C32 * j12 + T.C33 * j13 + c43 * j14;
-
- // const fvec cj02 = T.C20*j22 + T.C30*j23 + T.C40*j24;
- // const fvec cj12 = T.C21*j22 + T.C31*j23 + T.C41*j24;
- const fvec cj22 = T.C22 * j22 + T.C32 * j23 + c42 * j24;
- const fvec cj32 = T.C32 * j22 + T.C33 * j23 + c43 * j24;
-
- // const fvec cj03 = T.C20*j32 + T.C30*j33 + T.C40*j34;
- // const fvec cj13 = T.C21*j32 + T.C31*j33 + T.C41*j34;
- const fvec cj23 = T.C22 * j32 + T.C32 * j33 + c42 * j34;
- const fvec cj33 = T.C32 * j32 + T.C33 * j33 + c43 * j34;
-
- T.C40 += c42 * j02 + c43 * j03 + T.C44 * j04; // cj40
- T.C41 += c42 * j12 + c43 * j13 + T.C44 * j14; // cj41
- T.C42 = c42 * j22 + c43 * j23 + T.C44 * j24; // cj42
- T.C43 = c42 * j32 + c43 * j33 + T.C44 * j34; // cj43
-
- T.C00 = cj00 + j02 * cj20 + j03 * cj30 + j04 * T.C40;
- T.C10 = cj01 + j02 * cj21 + j03 * cj31 + j04 * T.C41;
- T.C11 = cj11 + j12 * cj21 + j13 * cj31 + j14 * T.C41;
-
- T.C20 = j22 * cj20 + j23 * cj30 + j24 * T.C40;
- T.C30 = j32 * cj20 + j33 * cj30 + j34 * T.C40;
- T.C21 = j22 * cj21 + j23 * cj31 + j24 * T.C41;
- T.C31 = j32 * cj21 + j33 * cj31 + j34 * T.C41;
- T.C22 = j22 * cj22 + j23 * cj32 + j24 * T.C42;
- T.C32 = j32 * cj22 + j33 * cj32 + j34 * T.C42;
- T.C33 = j32 * cj23 + j33 * cj33 + j34 * T.C43;
- //cout<<"Extrapolation ok"<<endl;
-}
-
-
-void L1ExtrapolateRungeKutta // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
- (L1TrackPar& T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
- fvec z_out, // extrapolate to this z position
- fvec qp0, // use Q/p linearisation at this value
- const L1FieldRegion& F, const fvec* w)
-{
- //
- // Forth-order Runge-Kutta method for solution of the equation
- // of motion of a particle with parameter qp = Q /P
- // in the magnetic field B()
- //
- // | x | tx
- // | y | ty
- // d/dz = | tx| = ft * ( ty * ( B(3)+tx*b(1) ) - (1+tx**2)*B(2) )
- // | ty| ft * (-tx * ( B(3)+ty+b(2) - (1+ty**2)*B(1) ) ,
- //
- // where ft = c_light*qp*sqrt ( 1 + tx**2 + ty**2 ) .
- //
- // In the following for RK step :
- //
- // x=x[0], y=x[1], tx=x[3], ty=x[4].
- //
- //========================================================================
- //
- // NIM A395 (1997) 169-184; NIM A426 (1999) 268-282.
- //
- // the routine is based on LHC(b) utility code
- //
- //========================================================================
-
-
- cnst c_light(0.000299792458);
-
- const fvec a[4] = {fvec(0.), fvec(0.5), fvec(0.5), fvec(1.)};
- const fvec c[4] = {fvec(1. / 6.), fvec(1. / 3.), fvec(1. / 3.), fvec(1. / 6.)};
- const fvec b[4] = {fvec(0.), fvec(0.5), fvec(0.5), fvec(1.)};
-
- int step4;
- fvec k[16], x0[4], x[4], k1[16];
- fvec Ax[4], Ay[4], Ax_tx[4], Ay_tx[4], Ax_ty[4], Ay_ty[4];
-
- //----------------------------------------------------------------
-
- if (w) { z_out = iif(fvec::Zero() < *w, z_out, T.z); }
-
- fvec qp_in = T.qp;
- const fvec z_in = T.z;
- const fvec h = (z_out - T.z);
- fvec hC = h * c_light;
- x0[0] = T.x;
- x0[1] = T.y;
- x0[2] = T.tx;
- x0[3] = T.ty;
- //
- // Runge-Kutta step
- //
-
- int step;
- int i;
-
- for (step = 0; step < 4; ++step) {
-
- for (i = 0; i < 4; ++i) {
- if (step == 0) { x[i] = x0[i]; }
- else {
- x[i] = x0[i] + b[step] * k[step * 4 - 4 + i];
- }
- }
- fvec B[3];
- F.Get(x[0], x[1], z_in + a[step] * h, B);
-
- fvec tx = x[2];
- fvec ty = x[3];
- fvec tx2 = tx * tx;
- fvec ty2 = ty * ty;
- fvec txty = tx * ty;
- fvec tx2ty21 = fvec(1.) + tx2 + ty2;
- // if( tx2ty21 > 1.e4 ) return 1;
- fvec I_tx2ty21 = fvec(1.) / tx2ty21 * qp0;
- fvec tx2ty2 = sqrt(tx2ty21);
- // fvec I_tx2ty2 = qp0 * hC / tx2ty2 ; unsused ???
- tx2ty2 *= hC;
- fvec tx2ty2qp = tx2ty2 * qp0;
- Ax[step] = (txty * B[0] + ty * B[2] - (fvec(1.) + tx2) * B[1]) * tx2ty2;
- Ay[step] = (-txty * B[1] - tx * B[2] + (fvec(1.) + ty2) * B[0]) * tx2ty2;
-
- Ax_tx[step] = Ax[step] * tx * I_tx2ty21 + (ty * B[0] - fvec(2.) * tx * B[1]) * tx2ty2qp;
- Ax_ty[step] = Ax[step] * ty * I_tx2ty21 + (tx * B[0] + B[2]) * tx2ty2qp;
- Ay_tx[step] = Ay[step] * tx * I_tx2ty21 + (-ty * B[1] - B[2]) * tx2ty2qp;
- Ay_ty[step] = Ay[step] * ty * I_tx2ty21 + (-tx * B[1] + fvec(2.) * ty * B[0]) * tx2ty2qp;
-
- step4 = step * 4;
- k[step4] = tx * h;
- k[step4 + 1] = ty * h;
- k[step4 + 2] = Ax[step] * qp0;
- k[step4 + 3] = Ay[step] * qp0;
-
- } // end of Runge-Kutta steps
-
- {
- T.x = x0[0] + c[0] * k[0] + c[1] * k[4 + 0] + c[2] * k[8 + 0] + c[3] * k[12 + 0];
- T.y = x0[1] + c[0] * k[1] + c[1] * k[4 + 1] + c[2] * k[8 + 1] + c[3] * k[12 + 1];
- T.tx = x0[2] + c[0] * k[2] + c[1] * k[4 + 2] + c[2] * k[8 + 2] + c[3] * k[12 + 2];
- T.ty = x0[3] + c[0] * k[3] + c[1] * k[4 + 3] + c[2] * k[8 + 3] + c[3] * k[12 + 3];
- T.z = z_out;
- }
-
- //
- // Derivatives dx/dqp
- //
-
- x0[0] = fvec(0.);
- x0[1] = fvec(0.);
- x0[2] = fvec(0.);
- x0[3] = fvec(0.);
-
- //
- // Runge-Kutta step for derivatives dx/dqp
-
- for (step = 0; step < 4; ++step) {
- for (i = 0; i < 4; ++i) {
- if (step == 0) { x[i] = x0[i]; }
- else {
- x[i] = x0[i] + b[step] * k1[step * 4 - 4 + i];
- }
- }
- step4 = step * 4;
- k1[step4] = x[2] * h;
- k1[step4 + 1] = x[3] * h;
- k1[step4 + 2] = Ax[step] + Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
- k1[step4 + 3] = Ay[step] + Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
-
- } // end of Runge-Kutta steps for derivatives dx/dqp
-
- fvec J[25];
-
- for (i = 0; i < 4; ++i) {
- J[20 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[4 + i] + c[2] * k1[8 + i] + c[3] * k1[12 + i];
- }
- J[24] = fvec(1.);
- //
- // end of derivatives dx/dqp
- //
-
- // Derivatives dx/tx
- //
-
- x0[0] = fvec(0.);
- x0[1] = fvec(0.);
- x0[2] = fvec(1.);
- x0[3] = fvec(0.);
-
- //
- // Runge-Kutta step for derivatives dx/dtx
- //
-
- for (step = 0; step < 4; ++step) {
- for (i = 0; i < 4; ++i) {
- if (step == 0) { x[i] = x0[i]; }
- else if (i != 2) {
- x[i] = x0[i] + b[step] * k1[step * 4 - 4 + i];
- }
- }
- step4 = step * 4;
- k1[step4] = x[2] * h;
- k1[step4 + 1] = x[3] * h;
- // k1[step4+2] = Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
- k1[step4 + 3] = Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
-
- } // end of Runge-Kutta steps for derivatives dx/dtx
-
- for (i = 0; i < 4; ++i) {
- if (i != 2) { J[10 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[4 + i] + c[2] * k1[8 + i] + c[3] * k1[12 + i]; }
- }
- // end of derivatives dx/dtx
- J[12] = fvec(1.);
- J[14] = fvec(0.);
-
- // Derivatives dx/ty
- //
-
- x0[0] = fvec(0.);
- x0[1] = fvec(0.);
- x0[2] = fvec(0.);
- x0[3] = fvec(1.);
-
- //
- // Runge-Kutta step for derivatives dx/dty
- //
-
- for (step = 0; step < 4; ++step) {
- for (i = 0; i < 4; ++i) {
- if (step == 0) {
- x[i] = x0[i]; // ty fixed
- }
- else if (i != 3) {
- x[i] = x0[i] + b[step] * k1[step * 4 - 4 + i];
- }
- }
- step4 = step * 4;
- k1[step4] = x[2] * h;
- k1[step4 + 1] = x[3] * h;
- k1[step4 + 2] = Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
- // k1[step4+3] = Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
-
- } // end of Runge-Kutta steps for derivatives dx/dty
-
- for (i = 0; i < 3; ++i) {
- J[15 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[4 + i] + c[2] * k1[8 + i] + c[3] * k1[12 + i];
- }
- // end of derivatives dx/dty
- J[18] = fvec(1.);
- J[19] = fvec(0.);
-
- //
- // derivatives dx/dx and dx/dy
-
- for (i = 0; i < 10; ++i) {
- J[i] = fvec(0.);
- }
- J[0] = fvec(1.);
- J[6] = fvec(1.);
-
- fvec dqp = qp_in - qp0;
-
- { // update parameters
- T.x += J[5 * 4 + 0] * dqp;
- T.y += J[5 * 4 + 1] * dqp;
- T.tx += J[5 * 4 + 2] * dqp;
- T.ty += J[5 * 4 + 3] * dqp;
- }
-
- // covariance matrix transport
-
- // // if ( C_in&&C_out ) CbmKFMath::multQtSQ( 5, J, C_in, C_out); // TODO
- // j(0,2) = J[5*2 + 0];
- // j(1,2) = J[5*2 + 1];
- // j(2,2) = J[5*2 + 2];
- // j(3,2) = J[5*2 + 3];
- //
- // j(0,3) = J[5*3 + 0];
- // j(1,3) = J[5*3 + 1];
- // j(2,3) = J[5*3 + 2];
- // j(3,3) = J[5*3 + 3];
- //
- // j(0,4) = J[5*4 + 0];
- // j(1,4) = J[5*4 + 1];
- // j(2,4) = J[5*4 + 2];
- // j(3,4) = J[5*4 + 3];
-
- const fvec c42 = T.C42, c43 = T.C43;
-
- const fvec cj00 = T.C00 + T.C20 * J[5 * 2 + 0] + T.C30 * J[5 * 3 + 0] + T.C40 * J[5 * 4 + 0];
- // const fvec cj10 = T.C10 + T.C21*J[5*2 + 0] + T.C31*J[5*3 + 0] + T.C41*J[5*4 + 0];
- const fvec cj20 = T.C20 + T.C22 * J[5 * 2 + 0] + T.C32 * J[5 * 3 + 0] + c42 * J[5 * 4 + 0];
- const fvec cj30 = T.C30 + T.C32 * J[5 * 2 + 0] + T.C33 * J[5 * 3 + 0] + c43 * J[5 * 4 + 0];
-
- const fvec cj01 = T.C10 + T.C20 * J[5 * 2 + 1] + T.C30 * J[5 * 3 + 1] + T.C40 * J[5 * 4 + 1];
- const fvec cj11 = T.C11 + T.C21 * J[5 * 2 + 1] + T.C31 * J[5 * 3 + 1] + T.C41 * J[5 * 4 + 1];
- const fvec cj21 = T.C21 + T.C22 * J[5 * 2 + 1] + T.C32 * J[5 * 3 + 1] + c42 * J[5 * 4 + 1];
- const fvec cj31 = T.C31 + T.C32 * J[5 * 2 + 1] + T.C33 * J[5 * 3 + 1] + c43 * J[5 * 4 + 1];
-
- // const fvec cj02 = T.C20*J[5*2 + 2] + T.C30*J[5*3 + 2] + T.C40*J[5*4 + 2];
- // const fvec cj12 = T.C21*J[5*2 + 2] + T.C31*J[5*3 + 2] + T.C41*J[5*4 + 2];
- const fvec cj22 = T.C22 * J[5 * 2 + 2] + T.C32 * J[5 * 3 + 2] + c42 * J[5 * 4 + 2];
- const fvec cj32 = T.C32 * J[5 * 2 + 2] + T.C33 * J[5 * 3 + 2] + c43 * J[5 * 4 + 2];
-
- // const fvec cj03 = T.C20*J[5*2 + 3] + T.C30*J[5*3 + 3] + T.C40*J[5*4 + 3];
- // const fvec cj13 = T.C21*J[5*2 + 3] + T.C31*J[5*3 + 3] + T.C41*J[5*4 + 3];
- const fvec cj23 = T.C22 * J[5 * 2 + 3] + T.C32 * J[5 * 3 + 3] + c42 * J[5 * 4 + 3];
- const fvec cj33 = T.C32 * J[5 * 2 + 3] + T.C33 * J[5 * 3 + 3] + c43 * J[5 * 4 + 3];
-
- T.C40 += c42 * J[5 * 2 + 0] + c43 * J[5 * 3 + 0] + T.C44 * J[5 * 4 + 0]; // cj40
- T.C41 += c42 * J[5 * 2 + 1] + c43 * J[5 * 3 + 1] + T.C44 * J[5 * 4 + 1]; // cj41
- T.C42 = c42 * J[5 * 2 + 2] + c43 * J[5 * 3 + 2] + T.C44 * J[5 * 4 + 2];
- T.C43 = c42 * J[5 * 2 + 3] + c43 * J[5 * 3 + 3] + T.C44 * J[5 * 4 + 3];
-
- T.C00 = cj00 + J[5 * 2 + 0] * cj20 + J[5 * 3 + 0] * cj30 + J[5 * 4 + 0] * T.C40;
- T.C10 = cj01 + J[5 * 2 + 0] * cj21 + J[5 * 3 + 0] * cj31 + J[5 * 4 + 0] * T.C41;
- T.C11 = cj11 + J[5 * 2 + 1] * cj21 + J[5 * 3 + 1] * cj31 + J[5 * 4 + 1] * T.C41;
-
- T.C20 = J[5 * 2 + 2] * cj20 + J[5 * 3 + 2] * cj30 + J[5 * 4 + 2] * T.C40;
- T.C30 = J[5 * 2 + 3] * cj20 + J[5 * 3 + 3] * cj30 + J[5 * 4 + 3] * T.C40;
- T.C21 = J[5 * 2 + 2] * cj21 + J[5 * 3 + 2] * cj31 + J[5 * 4 + 2] * T.C41;
- T.C31 = J[5 * 2 + 3] * cj21 + J[5 * 3 + 3] * cj31 + J[5 * 4 + 3] * T.C41;
- T.C22 = J[5 * 2 + 2] * cj22 + J[5 * 3 + 2] * cj32 + J[5 * 4 + 2] * T.C42;
- T.C32 = J[5 * 2 + 3] * cj22 + J[5 * 3 + 3] * cj32 + J[5 * 4 + 3] * T.C42;
- T.C33 = J[5 * 2 + 3] * cj23 + J[5 * 3 + 3] * cj33 + J[5 * 4 + 3] * T.C43;
-}
-
-
-void L1Extrapolate0(L1TrackPar& T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
- fvec z_out, // extrapolate to this z position
- L1FieldRegion& F) // magneti
-{
- //
- // Part of the analytic extrapolation formula with error (c_light*B*dz)^4/4!
- //
-
- cnst c_light(0.000299792458);
-
- cnst c1(1.), c2i(1. / 2.), c3i(1. / 3.), c6i(1. / 6.), c12i(1. / 12.);
-
- const fvec qp = T.qp;
-
-
- const fvec dz = (z_out - T.z);
- const fvec dz2 = dz * dz;
-
- // construct coefficients
-
- const fvec x = T.tx;
- const fvec y = T.ty;
- // const fvec z = T.z;
- const fvec xx = x * x;
- const fvec xy = x * y;
- const fvec yy = y * y;
- const fvec Ay = -xx - c1;
- const fvec Bx = yy + c1;
- const fvec ct = c_light * sqrt(c1 + xx + yy);
-
- const fvec dzc2i = dz * c2i;
- const fvec dz2c3i = dz2 * c3i;
- const fvec dzc6i = dz * c6i;
- const fvec dz2c12i = dz2 * c12i;
-
- const fvec sx = (F.cx0 + F.cx1 * dzc2i + F.cx2 * dz2c3i);
- const fvec sy = (F.cy0 + F.cy1 * dzc2i + F.cy2 * dz2c3i);
- const fvec sz = (F.cz0 + F.cz1 * dzc2i + F.cz2 * dz2c3i);
-
- const fvec Sx = (F.cx0 * c2i + F.cx1 * dzc6i + F.cx2 * dz2c12i);
- const fvec Sy = (F.cy0 * c2i + F.cy1 * dzc6i + F.cy2 * dz2c12i);
- const fvec Sz = (F.cz0 * c2i + F.cz1 * dzc6i + F.cz2 * dz2c12i);
-
-
- const fvec ctdz = ct * dz;
- const fvec ctdz2 = ct * dz2;
-
- const fvec j04 = ctdz2 * (Sx * xy + Sy * Ay + Sz * y);
- const fvec j14 = ctdz2 * (Sx * Bx - Sy * xy - Sz * x);
- const fvec j24 = ctdz * (sx * xy + sy * Ay + sz * y);
- const fvec j34 = ctdz * (sx * Bx - sy * xy - sz * x);
-
- T.x += x * dz + j04 * qp;
- T.y += y * dz + j14 * qp;
- T.tx += j24 * qp;
- T.ty += j34 * qp;
- T.z += dz;
- //T.t += sqrt((x-T.x)*(x-T.x)+(y-T.y)*(y-T.y)+dz*dz)/29.9792458f;
-
- // covariance matrix transport
-
- const fvec cj00 = T.C00 + T.C20 * dz + T.C40 * j04;
- // const fvec cj10 = T.C10 + T.C21*dz + T.C41*j04;
- const fvec cj20 = T.C20 + T.C22 * dz + T.C42 * j04;
- const fvec cj30 = T.C30 + T.C32 * dz + T.C43 * j04;
-
- const fvec cj01 = T.C10 + T.C30 * dz + T.C40 * j14;
- const fvec cj11 = T.C11 + T.C31 * dz + T.C41 * j14;
- const fvec cj21 = T.C21 + T.C32 * dz + T.C42 * j14;
- const fvec cj31 = T.C31 + T.C33 * dz + T.C43 * j14;
-
- // const fvec cj02 = T.C20 + T.C40*j24;
- // const fvec cj12 = T.C21 + T.C41*j24;
- const fvec cj22 = T.C22 + T.C42 * j24;
- const fvec cj32 = T.C32 + T.C43 * j24;
-
- // const fvec cj03 = T.C30 + T.C40*j34;
- // const fvec cj13 = T.C31 + T.C41*j34;
- // const fvec cj23 = T.C32 + T.C42*j34;
- const fvec cj33 = T.C33 + T.C43 * j34;
-
- T.C40 += T.C42 * dz + T.C44 * j04; // cj40
- T.C41 += T.C43 * dz + T.C44 * j14; // cj41
- T.C42 += T.C44 * j24; // cj42
- T.C43 += T.C44 * j34; // cj43
-
- T.C00 = cj00 + dz * cj20 + j04 * T.C40;
- T.C10 = cj01 + dz * cj21 + j04 * T.C41;
- T.C11 = cj11 + dz * cj31 + j14 * T.C41;
-
- T.C20 = cj20 + j24 * T.C40;
- T.C30 = cj30 + j34 * T.C40;
- T.C21 = cj21 + j24 * T.C41;
- T.C31 = cj31 + j34 * T.C41;
- T.C22 = cj22 + j24 * T.C42;
- T.C32 = cj32 + j34 * T.C42;
- T.C33 = cj33 + j34 * T.C43;
-}
-
-
-void L1ExtrapolateTime(L1TrackPar& T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
- fvec dz, // extrapolate to this z position
- fvec timeInfo)
-{
- cnst c_light(29.9792458);
- dz.setZero(timeInfo <= fvec::Zero());
-
- T.t += dz * sqrt(T.tx * T.tx + T.ty * T.ty + fvec(1.)) / c_light;
-
- const fvec k1 = dz * T.tx / (c_light * sqrt(T.tx * T.tx + T.ty * T.ty + fvec(1.)));
- const fvec k2 = dz * T.ty / (c_light * sqrt(T.tx * T.tx + T.ty * T.ty + fvec(1.)));
-
- fvec c52 = T.C52;
- fvec c53 = T.C53;
-
- T.C50 += k1 * T.C20 + k2 * T.C30;
- T.C51 += k1 * T.C21 + k2 * T.C31;
- T.C52 += k1 * T.C22 + k2 * T.C32;
- T.C53 += k1 * T.C32 + k2 * T.C33;
- T.C54 += k1 * T.C42 + k2 * T.C43;
- T.C55 += k1 * (T.C52 + c52) + k2 * (T.C53 + c53);
-}
-
-void L1ExtrapolateLine(L1TrackPar& T, fvec z_out)
-{
- fvec dz = z_out - T.z;
-
- T.x += T.tx * dz;
- T.y += T.ty * dz;
- T.z += dz;
-
- const fvec dzC32_in = dz * T.C32;
-
- T.C21 += dzC32_in;
- T.C10 += dz * (T.C21 + T.C30);
-
- const fvec C20_in = T.C20;
-
- T.C20 += dz * T.C22;
- T.C00 += dz * (T.C20 + C20_in);
-
- const fvec C31_in = T.C31;
-
- T.C31 += dz * T.C33;
- T.C11 += dz * (T.C31 + C31_in);
- T.C30 += dzC32_in;
-
- T.C40 += dz * T.C42;
- T.C41 += dz * T.C43;
-}
-
-
-void L1ExtrapolateJXY // is not used currently
- (fvec& tx, fvec& ty,
- fvec& qp, // input track parameters
- fvec dz, // extrapolate to this dz
- L1FieldRegion& F, fvec& extrDx, fvec& extrDy, fvec extrJ[])
-{
- //
- // Part of the analytic extrapolation formula with error (c_light*B*dz)^4/4!
- //
-
- cnst c_light(0.000299792458);
-
- cnst c1(1.), c2(2.), c3(3.), c4(4.), c6(6.), c9(9.), c15(15.), c18(18.), c45(45.), c2i(1. / 2.), c6i(1. / 6.),
- c12i(1. / 12.);
-
- fvec dz2 = dz * dz;
- fvec dz3 = dz2 * dz;
-
- // construct coefficients
-
- fvec x = tx;
- fvec y = ty;
- fvec xx = x * x;
- fvec xy = x * y;
- fvec yy = y * y;
- fvec y2 = y * c2;
- fvec x2 = x * c2;
- fvec x4 = x * c4;
- fvec xx31 = xx * c3 + c1;
- fvec xx159 = xx * c15 + c9;
-
- fvec Ay = -xx - c1;
- fvec Ayy = x * (xx * c3 + c3);
- fvec Ayz = -c2 * xy;
- fvec Ayyy = -(c15 * xx * xx + c18 * xx + c3);
-
- fvec Ayy_x = c3 * xx31;
- fvec Ayyy_x = -x4 * xx159;
-
- fvec Bx = yy + c1;
- fvec Byy = y * xx31;
- fvec Byz = c2 * xx + c1;
- fvec Byyy = -xy * xx159;
-
- fvec Byy_x = c6 * xy;
- fvec Byyy_x = -y * (c45 * xx + c9);
- fvec Byyy_y = -x * xx159;
-
- // end of coefficients calculation
-
- fvec t2 = c1 + xx + yy;
- fvec t = sqrt(t2);
- fvec h = qp * c_light;
- fvec ht = h * t;
-
- // get field integrals
- fvec Fx0 = F.cx0;
- fvec Fx1 = F.cx1 * dz;
- fvec Fx2 = F.cx2 * dz2;
- fvec Fy0 = F.cy0;
- fvec Fy1 = F.cy1 * dz;
- fvec Fy2 = F.cy2 * dz2;
- fvec Fz0 = F.cz0;
- fvec Fz1 = F.cz1 * dz;
- fvec Fz2 = F.cz2 * dz2;
-
- //
-
- fvec Sx = (Fx0 * c2i + Fx1 * c6i + Fx2 * c12i);
- fvec Sy = (Fy0 * c2i + Fy1 * c6i + Fy2 * c12i);
- fvec Sz = (Fz0 * c2i + Fz1 * c6i + Fz2 * c12i);
-
- fvec Syz;
- {
- constexpr double d = 1. / 2520.;
- cnst c00(21. * 20. * d), c01(21. * 5. * d), c02(21. * 2. * d), c10(7. * 30. * d), c11(7. * 9. * d),
- c12(7. * 4. * d), c20(2. * 63. * d), c21(2. * 21. * d), c22(2. * 10. * d);
- Syz = Fy0 * (c00 * Fz0 + c01 * Fz1 + c02 * Fz2) + Fy1 * (c10 * Fz0 + c11 * Fz1 + c12 * Fz2)
- + Fy2 * (c20 * Fz0 + c21 * Fz1 + c22 * Fz2);
- }
-
- fvec Syy;
- {
- constexpr double d = 1. / 2520.;
- cnst c00(420. * d), c01(21. * 15. * d), c02(21. * 8. * d), c03(63. * d), c04(70. * d), c05(20. * d);
- Syy = Fy0 * (c00 * Fy0 + c01 * Fy1 + c02 * Fy2) + Fy1 * (c03 * Fy1 + c04 * Fy2) + c05 * Fy2 * Fy2;
- }
-
- fvec Syyy;
- {
- constexpr double d = 181440.;
- cnst c000(7560 / d), c001(9 * 1008 / d), c002(5 * 1008 / d), c011(21 * 180 / d), c012(24 * 180 / d),
- c022(7 * 180 / d), c111(540 / d), c112(945 / d), c122(560 / d), c222(112 / d);
- fvec Fy22 = Fy2 * Fy2;
- Syyy = Fy0 * (Fy0 * (c000 * Fy0 + c001 * Fy1 + c002 * Fy2) + Fy1 * (c011 * Fy1 + c012 * Fy2) + c022 * Fy22)
- + Fy1 * (Fy1 * (c111 * Fy1 + c112 * Fy2) + c122 * Fy22) + c222 * Fy22 * Fy2;
- }
-
- fvec SA1 = Sx * xy + Sy * Ay + Sz * y;
- fvec SA1_x = Sx * y - Sy * x2;
- fvec SA1_y = Sx * x + Sz;
-
- fvec SB1 = Sx * Bx - Sy * xy - Sz * x;
- fvec SB1_x = -Sy * y - Sz;
- fvec SB1_y = Sx * y2 - Sy * x;
-
- fvec SA2 = Syy * Ayy + Syz * Ayz;
- fvec SA2_x = Syy * Ayy_x - Syz * y2;
- fvec SA2_y = -Syz * x2;
- fvec SB2 = Syy * Byy + Syz * Byz;
- fvec SB2_x = Syy * Byy_x + Syz * x4;
- fvec SB2_y = Syy * xx31;
-
- fvec SA3 = Syyy * Ayyy;
- fvec SA3_x = Syyy * Ayyy_x;
- fvec SB3 = Syyy * Byyy;
- fvec SB3_x = Syyy * Byyy_x;
- fvec SB3_y = Syyy * Byyy_y;
-
- fvec ht1 = ht * dz;
- fvec ht2 = ht * ht * dz2;
- fvec ht3 = ht * ht * ht * dz3;
- fvec ht1SA1 = ht1 * SA1;
- fvec ht1SB1 = ht1 * SB1;
- fvec ht2SA2 = ht2 * SA2;
- fvec ht2SB2 = ht2 * SB2;
- fvec ht3SA3 = ht3 * SA3;
- fvec ht3SB3 = ht3 * SB3;
-
- extrDx = (x + ht1SA1 + ht2SA2 + ht3SA3) * dz;
- extrDy = (y + ht1SB1 + ht2SB2 + ht3SB3) * dz;
-
- fvec ctdz2 = c_light * t * dz2;
-
- fvec t2i = c1 / t2;
- fvec xt2i = x * t2i;
- fvec yt2i = y * t2i;
- fvec tmp0 = ht1SA1 + c2 * ht2SA2 + c3 * ht3SA3;
- fvec tmp1 = ht1SB1 + c2 * ht2SB2 + c3 * ht3SB3;
-
- extrJ[0] = dz * (c1 + xt2i * tmp0 + ht1 * SA1_x + ht2 * SA2_x + ht3 * SA3_x); // j02
- extrJ[1] = dz * (yt2i * tmp0 + ht1 * SA1_y + ht2 * SA2_y); // j03
- extrJ[2] = ctdz2 * (SA1 + c2 * ht1 * SA2 + c3 * ht2 * SA3); // j04
- extrJ[3] = dz * (xt2i * tmp1 + ht1 * SB1_x + ht2 * SB2_x + ht3 * SB3_x); // j12
- extrJ[4] = dz * (c1 + yt2i * tmp1 + ht1 * SB1_y + ht2 * SB2_y + ht3 * SB3_y); // j13
- extrJ[5] = ctdz2 * (SB1 + c2 * ht1 * SB2 + c3 * ht2 * SB3); // j14
-}
-
-void L1ExtrapolateJXY0(fvec& tx,
- fvec& ty, // input track slopes
- fvec dz, // extrapolate to this dz position
- L1FieldRegion& F, fvec& extrDx, fvec& extrDy, fvec& J04, fvec& J14)
-{
- cnst c_light(0.000299792458), c1(1.), c2i(0.5), c6i(1. / 6.), c12i(1. / 12.);
-
- fvec dz2 = dz * dz;
- fvec dzc6i = dz * c6i;
- fvec dz2c12i = dz2 * c12i;
-
- fvec xx = tx * tx;
- fvec yy = ty * ty;
- fvec xy = tx * ty;
-
- fvec Ay = -xx - c1;
- fvec Bx = yy + c1;
-
- fvec ctdz2 = c_light * sqrt(c1 + xx + yy) * dz2;
-
- fvec Sx = F.cx0 * c2i + F.cx1 * dzc6i + F.cx2 * dz2c12i;
- fvec Sy = F.cy0 * c2i + F.cy1 * dzc6i + F.cy2 * dz2c12i;
- fvec Sz = F.cz0 * c2i + F.cz1 * dzc6i + F.cz2 * dz2c12i;
-
- extrDx = tx * dz;
- extrDy = ty * dz;
- J04 = ctdz2 * (Sx * xy + Sy * Ay + Sz * ty);
- J14 = ctdz2 * (Sx * Bx - Sy * xy - Sz * tx);
-}
-
-#undef cnst
diff --git a/reco/L1/L1Algo/L1Extrapolation.h b/reco/L1/L1Algo/L1Extrapolation.h
deleted file mode 100644
index f4ec3a4b5c8d5d8c53afe178ee69f9dd2dcee3a0..0000000000000000000000000000000000000000
--- a/reco/L1/L1Algo/L1Extrapolation.h
+++ /dev/null
@@ -1,60 +0,0 @@
-/* Copyright (C) 2007-2017 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Maksym Zyzak, Sergey Gorbunov [committer], Igor Kulakov */
-
-#ifndef L1Extrapolation_h
-#define L1Extrapolation_h
-
-#include "L1Def.h"
-#include "L1Field.h"
-#include "L1TrackPar.h"
-
-void L1Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
- (L1TrackPar& T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
- fvec z_out, // extrapolate to this z position
- fvec qp0, // use Q/p linearisation at this value
- const L1FieldRegion& F, const fvec* w = nullptr);
-
-void L1Extrapolate0(L1TrackPar& T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
- fvec z_out, // extrapolate to this z position
- L1FieldRegion& F);
-
-void L1ExtrapolateTime(L1TrackPar& T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
- fvec dz, // extrapolate to this z position
- fvec timeInfo = fvec::One());
-
-void L1ExtrapolateLine(L1TrackPar& T, fvec z_out);
-
-void L1ExtrapolateJXY // is not used currently
- (fvec& tx, fvec& ty,
- fvec& qp, // input track parameters
- fvec dz, // extrapolate to this dz
- L1FieldRegion& F, fvec& extrDx, fvec& extrDy, fvec extrJ[]);
-
-void L1ExtrapolateJXY0(fvec& tx,
- fvec& ty, // input track slopes
- fvec dz, // extrapolate to this dz position
- L1FieldRegion& F, fvec& extrDx, fvec& extrDy, fvec& J04, fvec& J14);
-
-
-inline void L1ExtrapolateXC00Line(const L1TrackPar& T, fvec z_out, fvec& x, fvec& C00)
-{
- const fvec dz = (z_out - T.z);
- x = T.x + T.tx * dz;
- C00 = T.C00 + dz * (2 * T.C20 + dz * T.C22);
-}
-
-inline void L1ExtrapolateYC11Line(const L1TrackPar& T, fvec z_out, fvec& y, fvec& C11)
-{
- const fvec dz = (z_out - T.z);
- y = T.y + T.ty * dz;
- C11 = T.C11 + dz * (2 * T.C31 + dz * T.C33);
-}
-
-inline void L1ExtrapolateC10Line(const L1TrackPar& T, fvec z_out, fvec& C10)
-{
- const fvec dz = (z_out - T.z);
- C10 = T.C10 + dz * (T.C21 + T.C30 + dz * T.C32);
-}
-
-#endif
diff --git a/reco/L1/L1Algo/L1Fit.cxx b/reco/L1/L1Algo/L1Fit.cxx
deleted file mode 100644
index ffa2a3ab59b81a8fa0e77651db5fc1d1741ec865..0000000000000000000000000000000000000000
--- a/reco/L1/L1Algo/L1Fit.cxx
+++ /dev/null
@@ -1,12 +0,0 @@
-/* Copyright (C) 2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-/// \file L1Fit.cxx
-/// \brief Implementation of the L1Fit class
-/// \author Sergey Gorbunov <se.gorbunov@gsi.de>
-/// \date 15.03.2021
-
-#include "L1Fit.h"
-
-//ClassImp(L1Fit);
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
deleted file mode 100644
index b317ea963b45e2b8fe55726b2ff1564f33648ec4..0000000000000000000000000000000000000000
--- a/reco/L1/L1Algo/L1Fit.h
+++ /dev/null
@@ -1,88 +0,0 @@
-/* Copyright (C) 2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-/// \file L1Fit.h
-/// \brief Definition of the L1Fit class
-/// \author Sergey Gorbunov <se.gorbunov@gsi.de>
-/// \date 15.03.2021
-
-#ifndef L1Fit_H
-#define L1Fit_H
-
-#include "L1Constants.h"
-#include "L1Def.h"
-#include "L1Material.h"
-#include "L1TrackPar.h"
-
-///
-/// A collection of fit routines used by the CA tracker
-///
-class L1Fit {
-
-public:
- /// constructor
- L1Fit() {};
-
- /// destructor
- ~L1Fit() {};
-
- /// set particle mass for the fit
- void SetParticleMass(fvec mass)
- {
- fMass = mass;
- fMass2 = mass * mass;
- }
-
- /// get the particle mass
- fvec GetParticleMass() const { return fMass; }
-
- /// get the particle mass squared
- fvec GetParticleMass2() const { return fMass2; }
-
- /// === routines for material
-
- fvec BetheBlochIron(const float qp);
-
- fvec BetheBlochCarbon(const float qp);
-
- fvec BetheBlochAl(const float qp);
-
- static fvec ApproximateBetheBloch(const fvec& bg2);
-
- static fvec ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const fvec& kp1, const fvec& kp2, const fvec& kp3,
- const fvec& kp4);
-
- void EnergyLossCorrection(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w);
-
- template<int atomicZ>
- void EnergyLossCorrection(float atomicA, float rho, float radLen, L1TrackPar& T, const fvec& radThick, fvec& qp0,
- fvec direction, fvec w);
-
- void EnergyLossCorrectionAl(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w);
-
-
- void EnergyLossCorrectionCarbon(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w);
-
-
- void EnergyLossCorrectionIron(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w);
-
- void L1AddMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w);
-
- void L1AddThickMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream);
-
- void L1AddPipeMaterial(L1TrackPar& T, fvec qp0, fvec w);
-
- void L1AddTargetMaterial(L1TrackPar& T, fvec qp0, fvec w);
-
- const fvec PipeRadThick = 7.87e-3f; // 0.7 mm Aluminium // TODO: !
- const fvec TargetRadThick = 3.73e-2f * 2; // 250 mum Gold // TODO: !
-
-private:
- fvec fMass {L1Constants::phys::kMuonMass}; // muon mass by default
- fvec fMass2 {fMass * fMass};
-
- //ClassDefNV(L1Fit, 0)
-};
-
-#endif
diff --git a/reco/L1/L1Algo/L1FitMaterial.cxx b/reco/L1/L1Algo/L1FitMaterial.cxx
deleted file mode 100644
index 8491999c77e9171a5ef109435d9847a6d96d4c7c..0000000000000000000000000000000000000000
--- a/reco/L1/L1Algo/L1FitMaterial.cxx
+++ /dev/null
@@ -1,516 +0,0 @@
-/* Copyright (C) 2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#include "L1Def.h"
-#include "L1Fit.h"
-#include "L1Material.h"
-#include "L1TrackPar.h"
-
-//#define cnst static const fvec
-#define cnst const fvec
-
-//const fvec PipeRadThick = 7.87e-3f; // 0.7 mm Aluminium
-//const fvec TargetRadThick = 3.73e-2f * 2; // 250 mum Gold
-
-fvec L1Fit::BetheBlochIron(const float qp)
-{
-
- float K = 0.000307075; // [GeV*cm^2/g]
- float z = (qp > 0.) ? 1 : -1.;
- float Z = 26;
- float A = 55.845;
-
- float M = 0.10565f;
- float p = fabs(1. / qp); //GeV
- float E = sqrt(M * M + p * p);
- float beta = p / E;
- float betaSq = beta * beta;
- float gamma = E / M;
- float gammaSq = gamma * gamma;
-
- float I = 260 * 1e-9; // GeV
-
- float me = 0.000511; // GeV
- float ratio = me / M;
- float Tmax = (2 * me * betaSq * gammaSq) / (1 + 2 * gamma * ratio + ratio * ratio);
-
- // density correction
- float dc = 0.;
- if (p > 0.5) { // for particles above 1 Gev
- float rho = 7.87;
- float hwp = 28.816 * sqrt(rho * Z / A) * 1e-9; // GeV
- dc = log(hwp / I) + log(beta * gamma) - 0.5;
- }
-
- return K * z * z * (Z / A) * (1. / betaSq) * (0.5 * log(2 * me * betaSq * gammaSq * Tmax / (I * I)) - betaSq - dc);
-}
-
-fvec L1Fit::BetheBlochCarbon(const float qp)
-{
-
- constexpr float K = 0.000307075; // GeV * g^-1 * cm^2
- float z = (qp > 0.) ? 1 : -1.;
- constexpr float Z = 6;
- constexpr float A = 12.011;
-
- constexpr float M = 0.10565f;
- float p = fabs(1. / qp); //GeV
- float E = sqrt(M * M + p * p);
- float beta = p / E;
- float betaSq = beta * beta;
- float gamma = E / M;
- float gammaSq = gamma * gamma;
-
- const float I = 16 * std::pow(6, 0.9) * 1e-9; // GeV mean excitation energy in eV
-
- constexpr float me = 0.000511; // GeV
- constexpr float ratio = me / M;
- float Tmax = (2 * me * betaSq * gammaSq) / (1 + 2 * gamma * ratio + ratio * ratio);
-
- // density correction
- float dc = 0.;
- if (p > 0.5) { // for particles above 1 Gev
- constexpr float rho = 2.265;
- const float hwp = 28.816 * sqrt(rho * Z / A) * 1e-9; // GeV
- dc = log(hwp / I) + log(beta * gamma) - 0.5;
- }
-
- return K * z * z * (Z / A) * (1. / betaSq) * (0.5 * log(2 * me * betaSq * gammaSq * Tmax / (I * I)) - betaSq - dc);
-}
-
-fvec L1Fit::BetheBlochAl(const float qp)
-{
-
- constexpr float K = 0.000307075; // GeV * g^-1 * cm^2
- float z = (qp > 0.) ? 1 : -1.;
- constexpr float Z = 13;
- constexpr float A = 26.981;
-
- constexpr float M = 0.10565f;
- float p = fabs(1. / qp); //GeV
- float E = sqrt(M * M + p * p);
- float beta = p / E;
- float betaSq = beta * beta;
- float gamma = E / M;
- float gammaSq = gamma * gamma;
-
- const float I = 16 * std::pow(6, 0.9) * 1e-9; // GeV mean excitation energy in eV
-
- constexpr float me = 0.000511; // GeV
- constexpr float ratio = me / M;
- float Tmax = (2 * me * betaSq * gammaSq) / (1 + 2 * gamma * ratio + ratio * ratio);
-
- // density correction
- float dc = 0.;
- if (p > 0.5) { // for particles above 1 Gev
- constexpr float rho = 2.70;
- const float hwp = 28.816 * sqrt(rho * Z / A) * 1e-9; // GeV
- dc = log(hwp / I) + log(beta * gamma) - 0.5;
- }
-
- return K * z * z * (Z / A) * (1. / betaSq) * (0.5 * log(2 * me * betaSq * gammaSq * Tmax / (I * I)) - betaSq - dc);
-}
-
-fvec L1Fit::ApproximateBetheBloch(const fvec& bg2)
-{
- //
- // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
- //
- // bg2 - (beta*gamma)^2
- // kp0 - density [g/cm^3]
- // kp1 - density effect first junction point
- // kp2 - density effect second junction point
- // kp3 - mean excitation energy [GeV]
- // kp4 - mean Z/A
- //
- // The default values for the kp* parameters are for silicon.
- // The returned value is in [GeV/(g/cm^2)].
- //
-
- const fvec kp0 = 2.33f;
- const fvec kp1 = 0.20f;
- const fvec kp2 = 3.00f;
- const fvec kp3 = 173e-9f;
- const fvec kp4 = 0.49848f;
-
- constexpr float mK = 0.307075e-3f; // [GeV*cm^2/g]
- constexpr float _2me = 1.022e-3f; // [GeV/c^2]
- const fvec rho = kp0;
- const fvec x0 = kp1 * 2.303f;
- const fvec x1 = kp2 * 2.303f;
- const fvec mI = kp3;
- const fvec mZA = kp4;
- const fvec maxT = _2me * bg2; // neglecting the electron mass
-
- //*** Density effect
- fvec d2(0.f);
- const fvec x = 0.5f * log(bg2);
- const fvec lhwI = log(28.816f * 1e-9f * sqrt(rho * mZA) / mI);
-
- fmask init = x > x1;
- d2 = iif(init, lhwI + x - 0.5f, fvec::Zero());
- const fvec r = (x1 - x) / (x1 - x0);
- init = (x > x0) & (x1 > x);
- d2 = iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
-
- return mK * mZA * (fvec(1.f) + bg2) / bg2
- * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2) - d2);
-}
-
-fvec L1Fit::ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const fvec& kp1, const fvec& kp2, const fvec& kp3,
- const fvec& kp4)
-{
- //
- // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
- //
- // bg2 - (beta*gamma)^2
- // kp0 - density [g/cm^3]
- // kp1 - density effect first junction point
- // kp2 - density effect second junction point
- // kp3 - mean excitation energy [GeV]
- // kp4 - mean Z/A
- //
- // The default values for the kp* parameters are for silicon.
- // The returned value is in [GeV/(g/cm^2)].
- //
-
- // const fvec &kp0 = 2.33f;
- // const fvec &kp1 = 0.20f;
- // const fvec &kp2 = 3.00f;
- // const fvec &kp3 = 173e-9f;
- // const fvec &kp4 = 0.49848f;
-
- constexpr float mK = 0.307075e-3f; // [GeV*cm^2/g]
- constexpr float _2me = 1.022e-3f; // [GeV/c^2]
- const fvec& rho = kp0;
- const fvec x0 = kp1 * 2.303f;
- const fvec x1 = kp2 * 2.303f;
- const fvec& mI = kp3;
- const fvec& mZA = kp4;
- const fvec maxT = _2me * bg2; // neglecting the electron mass
-
- //*** Density effect
- fvec d2(0.f);
- const fvec x = 0.5f * log(bg2);
- const fvec lhwI = log(28.816f * 1e-9f * sqrt(rho * mZA) / mI);
-
- fmask init = x > x1;
- d2 = iif(init, lhwI + x - 0.5f, fvec::Zero());
- const fvec r = (x1 - x) / (x1 - x0);
- init = (x > x0) & (x1 > x);
- d2 = iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
-
- return mK * mZA * (fvec(1.f) + bg2) / bg2
- * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2) - d2);
-}
-
-
-void L1Fit::EnergyLossCorrection(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w)
-{
- const fvec p2 = 1.f / (T.qp * T.qp);
- const fvec E2 = fMass2 + p2;
-
- // fvec qp =T.qp;
-
- const fvec bethe = ApproximateBetheBloch(p2 / fMass2);
-
- // fvec bethe2 = BetheBlochAl(qp[0]);
-
-
- fvec tr = sqrt(fvec(1.f) + T.tx * T.tx + T.ty * T.ty);
-
- fvec dE = bethe * radThick * tr * 9.34961f * 2.33f;
- // fvec dE2 = bethe2 * radThick*tr*2.265f* 2.70f;
-
- // dE = fabs(dE_);//dE2;//bethe * radThick*tr * 9.34961f*2.33f ;
-
-
- const fvec E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - fMass2));
- fmask ok = (corr == corr) & (fvec::Zero() < w);
- corr = iif(ok, corr, fvec::One());
-
- qp0 *= corr;
- // fvec dqp = CalcQpAfterEloss(qp[0], (direction*dE)[0], fMass2[0]);
- // qp0 = dqp;
- // T.qp = dqp;
- T.qp *= corr;
- T.C40 *= corr;
- T.C41 *= corr;
- T.C42 *= corr;
- T.C43 *= corr;
- T.C44 *= corr * corr;
-}
-
-template<int atomicZ>
-void L1Fit::EnergyLossCorrection(float atomicA, float rho, float radLen, L1TrackPar& T, const fvec& radThick, fvec& qp0,
- fvec direction, fvec w)
-{
- const fvec p2 = 1.f / (T.qp * T.qp);
- const fvec E2 = fMass2 + p2;
-
- fvec qp = T.qp;
-
- fvec i;
- if (atomicZ < 13) i = (12. * atomicZ + 7.) * 1.e-9;
- else
- i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
-
- const fvec bethe = ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
-
- // const fvec bethe = ApproximateBetheBloch( p2/fMass2 );
-
- // fvec bethe2 = BetheBlochAl(qp[0]);
-
- fvec tr = sqrt(fvec(1.f) + T.tx * T.tx + T.ty * T.ty);
-
- // fvec dE2 = 2*bethe * radThick*tr * 9.34961f*2.33f ;
- fvec dE = bethe * radThick * tr * radLen * rho;
-
- const fvec E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - fMass2));
- fmask ok = (corr == corr) & (fvec::Zero() < w);
- corr = iif(ok, corr, fvec::One());
-
- qp0 *= corr;
- // fvec dqp = CalcQpAfterEloss(qp[0], (direction*dE)[0], fMass2[0]);
- // qp0 = dqp;
- // T.qp = dqp;
- T.qp *= corr;
-
- fvec P = fvec(fabs(1. / qp[0])); // GeV
-
- fvec Z(atomicZ);
- fvec A(atomicA);
- fvec RHO(rho);
-
- fvec STEP = radThick * tr * radLen;
- fvec EMASS(0.511 * 1e-3); // GeV
-
- fvec BETA = P / sqrt(E2);
- fvec GAMMA = sqrt(E2) / fMass;
-
- // Calculate xi factor (KeV).
- fvec XI = (fvec(153.5) * Z * STEP * RHO) / (A * BETA * BETA);
-
- // Maximum energy transfer to atomic electron (KeV).
- fvec ETA = BETA * GAMMA;
- fvec ETASQ = ETA * ETA;
- fvec RATIO = EMASS / fMass;
- fvec F1 = fvec(2.) * EMASS * ETASQ;
- fvec F2 = fvec(1.) + fvec(2.) * RATIO * GAMMA + RATIO * RATIO;
- fvec EMAX = fvec(1e6) * F1 / F2;
-
- fvec DEDX2 = XI * EMAX * (fvec(1.) - (BETA * BETA / fvec(2.))) * fvec(1e-12);
-
- fvec P2 = P * P;
- fvec SDEDX = ((E2) *DEDX2) / (P2 * P2 * P2);
-
-
- // T.C40 *= corr;
- // T.C41 *= corr;
- // T.C42 *= corr;
- // T.C43 *= corr;
- // T.C44 *= corr*corr;
- T.C44 += abs(SDEDX);
-
- // T.C40 *= corr;
- // T.C41 *= corr;
- // T.C42 *= corr;
- // T.C43 *= corr;
- // T.C44 *= corr*corr;
-}
-
-void L1Fit::EnergyLossCorrectionAl(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w)
-{
- constexpr int atomicZ = 13;
- constexpr float atomicA = 26.981f;
- constexpr float rho = 2.70f;
- constexpr float radLen = 2.265f;
- EnergyLossCorrection<atomicZ>(atomicA, rho, radLen, T, radThick, qp0, direction, w);
-}
-
-void L1Fit::EnergyLossCorrectionCarbon(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w)
-{
- constexpr int atomicZ = 6;
- constexpr float atomicA = 12.011f;
- constexpr float rho = 2.265;
- constexpr float radLen = 18.76f;
- EnergyLossCorrection<atomicZ>(atomicA, rho, radLen, T, radThick, qp0, direction, w);
-}
-
-void L1Fit::EnergyLossCorrectionIron(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w)
-{
- constexpr int atomicZ = 26;
- constexpr float atomicA = 55.845f;
- constexpr float rho = 7.87;
- constexpr float radLen = 1.758f;
- EnergyLossCorrection<atomicZ>(atomicA, rho, radLen, T, radThick, qp0, direction, w);
-}
-
-
-void L1Fit::L1AddMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w)
-{
- cnst ONE = 1.;
-
- fvec tx = T.tx;
- fvec ty = T.ty;
- fvec txtx = tx * tx;
- fvec tyty = ty * ty;
- fvec txtx1 = txtx + ONE;
- fvec h = txtx + tyty;
- fvec t = sqrt(txtx1 + tyty);
- fvec h2 = h * h;
- fvec qp0t = qp0 * t;
-
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
-
- fvec s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+fMass2*qp0*qp0t)*radThick*s0*s0 );
- fvec a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
-
- T.C22 += w * txtx1 * a;
- T.C32 += w * tx * ty * a;
- T.C33 += w * (ONE + tyty) * a;
-}
-
-
-// inline void L1Fit::L1AddThickMaterial( L1TrackPar &T, fvec radThick, fvec qp0, fvec thickness=0, fvec w = 1, fvec mass2 = 0.10565f*0.10565f, bool fDownstream = 1 )
-// {
-// cnst ONE = 1.;
-//
-// fvec tx = T.tx;
-// fvec ty = T.ty;
-// fvec txtx = tx*tx;
-// fvec tyty = ty*ty;
-// fvec txtx1 = txtx + ONE;
-// fvec h = txtx + tyty;
-// fvec t = sqrt(txtx1 + tyty);
-// fvec h2 = h*h;
-// fvec qp0t = qp0*t;
-//
-// cnst c1=0.0136f, c2=c1*0.038f, c3=c2*0.5f, c4=-c3/2.0f, c5=c3/3.0f, c6=-c3/4.0f;
-//
-// fvec s0 = (c1+c2*log(radThick) + c3*h + h2*(c4 + c5*h +c6*h2) )*qp0t;
-// //fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
-// fvec a = ( (t+mass2*qp0*qp0t)*radThick*s0*s0 );
-//
-// fvec D = (fDownstream) ? 1. : -1.;
-// fvec T23 = (thickness * thickness) / 3.0;
-// fvec T2 = thickness / 2.0;
-//
-//
-//
-//
-// T.C00 += w*txtx1*a * T23;
-// T.C10 += w*tx*ty*a * T23;
-// T.C20 += w*txtx1*a * D * T2;
-// T.C30 += w*tx*ty*a * D * T2;
-//
-// T.C11 += w*(ONE+tyty)*a * T23;
-// T.C21 += w*tx*ty*a * D * T2;
-// T.C31 += w*(ONE+tyty)*a * D * T2;
-//
-// T.C22 += w*txtx1*a;
-// T.C32 += w*tx*ty*a;
-// T.C33 += w*(ONE+tyty)*a;
-//
-// }
-
-void L1Fit::L1AddThickMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream)
-{
- cnst ONE = 1.;
-
- fvec tx = T.tx;
- fvec ty = T.ty;
- fvec txtx = tx * tx;
- fvec tyty = ty * ty;
- fvec txtx1 = txtx + ONE;
- fvec h = txtx + tyty;
- fvec t = sqrt(txtx1 + tyty);
- fvec h2 = h * h;
- fvec qp0t = qp0 * t;
-
- cnst c1(0.0136), c2 = c1 * fvec(0.038), c3 = c2 * fvec(0.5), c4 = -c3 / fvec(2.0), c5 = c3 / fvec(3.0),
- c6 = -c3 / fvec(4.0);
-
- fvec s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+fMass2*qp0*qp0t)*radThick*s0*s0 );
- fvec a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
-
- fvec D = (fDownstream) ? fvec(1.) : fvec(-1.);
- fvec T23 = (thickness * thickness) / fvec(3.0);
- fvec T2 = thickness / fvec(2.0);
-
-
- T.C00 += w * txtx1 * a * T23;
- T.C10 += w * tx * ty * a * T23;
- T.C20 += w * txtx1 * a * D * T2;
- T.C30 += w * tx * ty * a * D * T2;
-
- T.C11 += w * (ONE + tyty) * a * T23;
- T.C21 += w * tx * ty * a * D * T2;
- T.C31 += w * (ONE + tyty) * a * D * T2;
-
- T.C22 += w * txtx1 * a;
- T.C32 += w * tx * ty * a;
- T.C33 += w * (ONE + tyty) * a;
-}
-
-
-void L1Fit::L1AddPipeMaterial(L1TrackPar& T, fvec qp0, fvec w)
-{
- cnst ONE = 1.f;
-
- // static const fscal RadThick=0.0009f;//0.5/18.76;
- // static const fscal logRadThick=log(RadThick);
- //const fscal RadThick=0.0009f;//0.5/18.76;
-
- const fscal logRadThick = log(PipeRadThick[0]);
- fvec tx = T.tx;
- fvec ty = T.ty;
- fvec txtx = tx * tx;
- fvec tyty = ty * ty;
- fvec txtx1 = txtx + ONE;
- fvec h = txtx + tyty;
- fvec t = sqrt(txtx1 + tyty);
- fvec h2 = h * h;
- fvec qp0t = qp0 * t;
-
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
- fvec s0 = (c1 + c2 * fvec(logRadThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+fMass2*qp0*qp0t)*RadThick*s0*s0 );
- fvec a = ((t + fMass2 * qp0 * qp0t) * PipeRadThick * s0 * s0);
-
- T.C22 += w * txtx1 * a;
- T.C32 += w * tx * ty * a;
- T.C33 += w * (ONE + tyty) * a;
-}
-
-void L1Fit::L1AddTargetMaterial(L1TrackPar& T, fvec qp0, fvec w)
-{
- cnst ONE = 1.f;
-
- const fscal logRadThick = log(TargetRadThick[0]);
- fvec tx = T.tx;
- fvec ty = T.ty;
- fvec txtx = tx * tx;
- fvec tyty = ty * ty;
- fvec txtx1 = txtx + ONE;
- fvec h = txtx + tyty;
- fvec t = sqrt(txtx1 + tyty);
- fvec h2 = h * h;
- fvec qp0t = qp0 * t;
-
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
- fvec s0 = (c1 + c2 * fvec(logRadThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+fMass2*qp0*qp0t)*RadThick*s0*s0 );
- fvec a = ((t + fMass2 * qp0 * qp0t) * TargetRadThick * s0 * s0);
-
- T.C22 += w * txtx1 * a;
- T.C32 += w * tx * ty * a;
- T.C33 += w * (ONE + tyty) * a;
-}
-
-#undef cnst