From c2fd68bd7868665a5e047bcef1c36005e67be1a7 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Tue, 6 Sep 2022 01:36:42 +0000
Subject: [PATCH] L1: SIMD vectors: remove custom operators
---
reco/L1/L1Algo/L1CATrackFinder.cxx | 14 +++---
reco/L1/L1Algo/L1CloneMerger.cxx | 8 +--
reco/L1/L1Algo/L1Extrapolation.cxx | 43 +++++++++-------
reco/L1/L1Algo/L1Filtration.h | 8 +--
reco/L1/L1Algo/L1FitMaterial.cxx | 40 ++++++++-------
reco/L1/L1Algo/L1TrackExtender.cxx | 8 +--
reco/L1/L1Algo/L1TrackFitter.cxx | 14 +++---
reco/L1/L1Algo/L1TrackParFit.cxx | 49 ++++++++++---------
.../CbmL1RichENNRingFinderParallel.cxx | 14 +++---
reco/L1/ParticleFinder/CbmL1PFFitter.cxx | 16 +++---
reco/L1/vectors/L1vec.h | 2 +
reco/L1/vectors/L1vecPseudo.h | 17 +------
reco/L1/vectors/L1vecVc.h | 15 ------
13 files changed, 115 insertions(+), 133 deletions(-)
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 6fa5e9b0cc..54212c1dc7 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -224,14 +224,14 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
T.ty = ty;
T.t = time;
- T.qp = 0.;
- T.C20 = T.C21 = 0;
- T.C30 = T.C31 = T.C32 = 0;
- T.C40 = T.C41 = T.C42 = T.C43 = 0;
- T.C50 = T.C51 = T.C52 = T.C53 = T.C54 = 0;
- T.C22 = T.C33 = fMaxSlopePV * fMaxSlopePV / 9.;
+ T.qp = fvec(0.);
+ T.C20 = T.C21 = fvec(0.);
+ T.C30 = T.C31 = T.C32 = fvec(0.);
+ T.C40 = T.C41 = T.C42 = T.C43 = fvec(0.);
+ T.C50 = T.C51 = T.C52 = T.C53 = T.C54 = fvec(0.);
+ T.C22 = T.C33 = fMaxSlopePV * fMaxSlopePV / fvec(9.);
if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) T.C22 = T.C33 = 10;
- T.C44 = fMaxInvMom / 3. * fMaxInvMom / 3.;
+ T.C44 = fMaxInvMom / fvec(3.) * fMaxInvMom / fvec(3.);
T.C55 = timeEr * timeEr;
if (fParameters.DevIsFitSingletsFromTarget()) { // TODO: doesn't work. Why?
diff --git a/reco/L1/L1Algo/L1CloneMerger.cxx b/reco/L1/L1Algo/L1CloneMerger.cxx
index 3cd1acb2b6..7305303779 100644
--- a/reco/L1/L1Algo/L1CloneMerger.cxx
+++ b/reco/L1/L1Algo/L1CloneMerger.cxx
@@ -170,17 +170,17 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
stam = staf - 1;
fvec zm = frAlgo.GetParameters()->GetStation(stam).z;
- fvec xm = 0.5 * (Tf.x + Tf.tx * (zm - Tf.z) + Tb.x + Tb.tx * (zm - Tb.z));
- fvec ym = 0.5 * (Tf.y + Tf.ty * (zm - Tf.z) + Tb.y + Tb.ty * (zm - Tb.z));
+ fvec xm = fvec(0.5) * (Tf.x + Tf.tx * (zm - Tf.z) + Tb.x + Tb.tx * (zm - Tb.z));
+ fvec ym = fvec(0.5) * (Tf.y + Tf.ty * (zm - Tf.z) + Tb.y + Tb.ty * (zm - Tb.z));
frAlgo.GetParameters()->GetStation(stam).fieldSlice.GetFieldValue(xm, ym, fBm);
fld.Set(fBb, Tb.z, fBm, zm, fBf, Tf.z);
- fvec zMiddle = 0.5 * (Tb.z + Tf.z);
+ fvec zMiddle = fvec(0.5) * (Tb.z + Tf.z);
L1Extrapolate(Tf, zMiddle, Tf.qp, fld);
L1Extrapolate(Tb, zMiddle, Tb.qp, fld);
- fvec Chi2Tracks = 0.f;
+ fvec Chi2Tracks(0.);
FilterTracks(&(Tf.x), &(Tf.C00), &(Tb.x), &(Tb.C00), 0, 0, &Chi2Tracks);
if (Chi2Tracks[0] > 50) continue;
diff --git a/reco/L1/L1Algo/L1Extrapolation.cxx b/reco/L1/L1Algo/L1Extrapolation.cxx
index cb2bb78fdf..735528a41b 100644
--- a/reco/L1/L1Algo/L1Extrapolation.cxx
+++ b/reco/L1/L1Algo/L1Extrapolation.cxx
@@ -110,16 +110,18 @@ void L1ExtrapolateAnalytic(L1TrackPar& T, // input track parameters (x,y,tx,ty,
fvec syz;
{
- cnst d = 1. / 360., 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;
+ constexpr double d = 1. / 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;
{
- cnst d = 1. / 2520., 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;
+ 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);
}
@@ -129,15 +131,16 @@ void L1ExtrapolateAnalytic(L1TrackPar& T, // input track parameters (x,y,tx,ty,
fvec Syy;
{
- cnst d = 1. / 2520., 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;
+ 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;
{
- cnst d = 1. / 181440., 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;
+ constexpr double d = 1. / 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;
@@ -686,13 +689,13 @@ void L1ExtrapolateTime(L1TrackPar& T, // input track parameters (x,y,tx,ty,Q/p)
fvec dz, // extrapolate to this z position
fvec timeInfo)
{
- cnst c_light = 29.9792458;
+ cnst c_light(29.9792458);
dz.setZero(timeInfo <= fvec::Zero());
T.t += dz * sqrt((T.tx * T.tx) + (T.ty * T.ty) + 1) / c_light;
- const fvec k1 = dz * T.tx / (c_light * sqrt((T.tx * T.tx) + (T.ty * T.ty) + 1.));
- const fvec k2 = dz * T.ty / (c_light * sqrt((T.tx * T.tx) + (T.ty * T.ty) + 1.));
+ 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;
@@ -809,23 +812,25 @@ void L1ExtrapolateJXY // is not used currently
fvec Syz;
{
- cnst d = 1. / 2520., 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;
+ 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;
{
- cnst d = 1. / 2520., 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;
+ 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;
{
- cnst d = 1. / 181440., 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;
+ constexpr double d = 1. / 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;
diff --git a/reco/L1/L1Algo/L1Filtration.h b/reco/L1/L1Algo/L1Filtration.h
index 1e28b13418..f9975ef476 100644
--- a/reco/L1/L1Algo/L1Filtration.h
+++ b/reco/L1/L1Algo/L1Filtration.h
@@ -256,7 +256,7 @@ inline void L1FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fve
HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
- wi = 1. / (info.sigma2 + HCH);
+ wi = fvec(1.) / (info.sigma2 + HCH);
zetawi = zeta * wi;
chi2 += zeta * zetawi;
@@ -273,7 +273,7 @@ inline void L1FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fve
inline void L1FilterVtx(L1TrackPar& T, fvec x, fvec y, const L1XYMeasurementInfo& info, fvec extrDx, fvec extrDy,
fvec J[])
{
- cnst TWO = 2.;
+ cnst TWO(2.);
fvec zeta0, zeta1, S00, S10, S11, si;
fvec F00, F10, F20, F30, F40, F01, F11, F21, F31, F41;
@@ -304,14 +304,14 @@ inline void L1FilterVtx(L1TrackPar& T, fvec x, fvec y, const L1XYMeasurementInfo
S10 = info.C10 + F10 + J[3] * F20 + J[4] * F30 + J[5] * F40;
S11 = info.C11 + F11 + J[3] * F21 + J[4] * F31 + J[5] * F41;
- si = 1. / (S00 * S11 - S10 * S10);
+ si = fvec(1.) / (S00 * S11 - S10 * S10);
//si = fvec(rcp(fvec((S00*S11 - S10*S10)[0])));
fvec S00tmp = S00;
S00 = si * S11;
S10 = -si * S10;
S11 = si * S00tmp;
- T.chi2 += zeta0 * zeta0 * S00 + 2. * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11;
+ T.chi2 += zeta0 * zeta0 * S00 + fvec(2.) * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11;
T.NDF += TWO;
K00 = F00 * S00 + F01 * S10;
diff --git a/reco/L1/L1Algo/L1FitMaterial.cxx b/reco/L1/L1Algo/L1FitMaterial.cxx
index a5b4586b3e..2804436807 100644
--- a/reco/L1/L1Algo/L1FitMaterial.cxx
+++ b/reco/L1/L1Algo/L1FitMaterial.cxx
@@ -279,32 +279,32 @@ void L1Fit::EnergyLossCorrection(float atomicA, float rho, float radLen, L1Track
// T.qp = dqp;
T.qp *= corr;
- float P = fabs(1. / qp[0]); // GeV
+ fvec P = fvec(fabs(1. / qp[0])); // GeV
- float Z = atomicZ;
- float A = atomicA;
- float RHO = rho;
+ fvec Z(atomicZ);
+ fvec A(atomicA);
+ fvec RHO(rho);
- fvec STEP = radThick * tr * radLen;
- fvec EMASS = 0.511 * 1e-3; // GeV
+ 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 = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
+ 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 = 2. * EMASS * ETASQ;
- fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
- fvec EMAX = 1e6 * F1 / F2;
+ 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 * (1. - (BETA * BETA / 2.)) * 1e-12;
+ fvec DEDX2 = XI * EMAX * (fvec(1.) - (BETA * BETA / fvec(2.))) * fvec(1e-12);
- float P2 = P * P;
+ fvec P2 = P * P;
fvec SDEDX = ((E2) *DEDX2) / (P2 * P2 * P2);
@@ -456,15 +456,16 @@ void L1Fit::L1AddThickMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w, f
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;
+ 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) ? 1. : -1.;
- fvec T23 = (thickness * thickness) / 3.0;
- fvec T2 = thickness / 2.0;
+ 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;
@@ -495,11 +496,12 @@ void L1Fit::L1AddHalfMaterial(L1TrackPar& T, const L1MaterialInfo& info, fvec qp
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;
+ 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 * (info.logRadThick + log(0.5)) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
+ fvec s0 = (c1 + c2 * (info.logRadThick + fvec(log(0.5))) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
//fvec a = ( (ONE+fMass2*qp0*qp0t)*info.RadThick*0.5*s0*s0 );
- fvec a = ((t + fMass2 * qp0 * qp0t) * info.RadThick * 0.5 * s0 * s0);
+ fvec a = ((t + fMass2 * qp0 * qp0t) * info.RadThick * fvec(0.5) * s0 * s0);
T.C22 += txtx1 * a;
T.C32 += tx * ty * a;
diff --git a/reco/L1/L1Algo/L1TrackExtender.cxx b/reco/L1/L1Algo/L1TrackExtender.cxx
index ae8298e166..bfd01b3a60 100644
--- a/reco/L1/L1Algo/L1TrackExtender.cxx
+++ b/reco/L1/L1Algo/L1TrackExtender.cxx
@@ -67,9 +67,9 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
auto [x2, y2] = sta2.ConvUVtoXY<fvec>(u2, v2);
// fvec z2 = hit2.z;
- fvec dzi = 1. / (z1 - z0);
+ fvec dzi = fvec(1.) / (z1 - z0);
- const fvec vINF = .1;
+ const fvec vINF = fvec(.1);
T.x = x0;
T.y = y0;
if (initParams) {
@@ -82,8 +82,8 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
T.t = hit0.t;
// T.t[0]=(hit0.t+hit1.t+hit2.t)/3;
- T.chi2 = 0.;
- T.NDF = 2.;
+ T.chi2 = fvec(0.);
+ T.NDF = fvec(2.);
T.C00 = sta0.XYInfo.C00;
T.C10 = sta0.XYInfo.C10;
T.C11 = sta0.XYInfo.C11;
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 6aeb45a1f2..0d63bb0736 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -19,10 +19,10 @@ using std::vector;
namespace NS_L1TrackFitter
{
- const fvec c_light = 0.000299792458, c_light_i = 1. / c_light;
- const fvec ZERO = 0.;
- const fvec ONE = 1.;
- const fvec vINF = 0.1;
+ const fvec c_light(0.000299792458), c_light_i(fvec(1.) / c_light);
+ const fvec ZERO(0.);
+ const fvec ONE(1.);
+ const fvec vINF(0.1);
} // namespace NS_L1TrackFitter
using namespace NS_L1TrackFitter;
@@ -81,7 +81,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
auto [x2, y2] = sta1.ConvUVtoXY<fvec>(u2, v2);
// fvec z2 = hit2.z;
- fvec dzi = 1. / (z1 - z0);
+ fvec dzi = fvec(1.) / (z1 - z0);
// const fvec vINF = .1;
T.x = x0;
@@ -1470,7 +1470,7 @@ void L1Algo::GuessVec(L1TrackPar& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fve
fvec L, L1;
t.x = (Ai0 * a0 + Ai1 * a1 + Ai3 * a2) * det;
t.tx = (Ai1 * a0 + Ai2 * a1 + Ai4 * a2) * det;
- fvec txtx1 = 1. + t.tx * t.tx;
+ fvec txtx1 = fvec(1.) + t.tx * t.tx;
L = (Ai3 * a0 + Ai4 * a1 + Ai5 * a2) * det / txtx1;
L1 = L * t.tx;
A1 = A1 + A3 * L1;
@@ -1545,7 +1545,7 @@ void L1Algo::GuessVec(L1TrackParFit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy,
fvec L, L1;
t.fx = (Ai0 * a0 + Ai1 * a1 + Ai3 * a2) * det;
t.ftx = (Ai1 * a0 + Ai2 * a1 + Ai4 * a2) * det;
- fvec txtx1 = 1. + t.ftx * t.ftx;
+ fvec txtx1 = fvec(1.) + t.ftx * t.ftx;
L = (Ai3 * a0 + Ai4 * a1 + Ai5 * a2) * det / txtx1;
L1 = L * t.ftx;
A1 = A1 + A3 * L1;
diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx
index 9890852328..7cd752c3a2 100644
--- a/reco/L1/L1Algo/L1TrackParFit.cxx
+++ b/reco/L1/L1Algo/L1TrackParFit.cxx
@@ -225,10 +225,10 @@ void L1TrackParFit::ExtrapolateLine(fvec z_out, fvec* w)
fx += dz * ftx;
fy += dz * fty;
fz += dz;
- ft += dz * sqrt(1. + ftx * ftx + fty * fty) / c_light;
+ ft += dz * sqrt(fvec(1.) + ftx * ftx + fty * fty) / c_light;
- const fvec k1 = dz * ftx / (c_light * sqrt(ftx * ftx + fty * fty + 1.));
- const fvec k2 = dz * fty / (c_light * sqrt(ftx * ftx + fty * fty + 1.));
+ const fvec k1 = dz * ftx / (c_light * sqrt(ftx * ftx + fty * fty + fvec(1.)));
+ const fvec k2 = dz * fty / (c_light * sqrt(ftx * ftx + fty * fty + fvec(1.)));
const fvec dzC32_in = dz * C32;
@@ -263,7 +263,7 @@ void L1TrackParFit::ExtrapolateLine(fvec z_out, fvec* w)
void L1TrackParFit::ExtrapolateLine1(fvec z_out, fvec* w, fvec v)
{
- cnst c_light = 29.9792458;
+ cnst c_light(29.9792458);
fvec dz = (z_out - fz);
if (w) { dz.setZero(*w <= fvec::Zero()); }
@@ -272,10 +272,10 @@ void L1TrackParFit::ExtrapolateLine1(fvec z_out, fvec* w, fvec v)
fy += dz * fty;
fz += dz;
- ft += dz * sqrt(1. + ftx * ftx + fty * fty) / (v * c_light);
+ ft += dz * sqrt(fvec(1.) + ftx * ftx + fty * fty) / (v * c_light);
- const fvec k1 = dz * ftx / ((v * c_light) * sqrt(ftx * ftx + fty * fty + 1.));
- const fvec k2 = dz * fty / ((v * c_light) * sqrt(ftx * ftx + fty * fty + 1.));
+ const fvec k1 = dz * ftx / ((v * c_light) * sqrt(ftx * ftx + fty * fty + fvec(1.)));
+ const fvec k2 = dz * fty / ((v * c_light) * sqrt(ftx * ftx + fty * fty + fvec(1.)));
const fvec dzC32_in = dz * C32;
@@ -741,15 +741,16 @@ void L1TrackParFit::L1AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thi
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;
+ 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+mass2*qp0*qp0t)*radThick*s0*s0 );
fvec a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
- fvec D = (fDownstream) ? 1. : -1.;
- fvec T23 = (thickness * thickness) / 3.0;
- fvec T2 = thickness / 2.0;
+ fvec D = (fDownstream) ? fvec(1.) : fvec(-1.);
+ fvec T23 = (thickness * thickness) / fvec(3.0);
+ fvec T2 = thickness / fvec(2.0);
C00 += w * txtx1 * a * T23;
C10 += w * tx * ty * a * T23;
@@ -845,33 +846,33 @@ void L1TrackParFit::EnergyLossCorrection(float atomicA, float rho, float radLen,
qp0 *= corr;
fqp *= corr;
- float P = fabs(1. / qp0[0]); // GeV
+ fvec P(fabs(1. / qp0[0])); // GeV
- float Z = atomicZ;
- float A = atomicA;
- float RHO = rho;
+ fvec Z(atomicZ);
+ fvec A(atomicA);
+ fvec RHO(rho);
- fvec STEP = radThick * tr * radLen;
- fvec EMASS = 0.511 * 1e-3; // GeV
+ fvec STEP = radThick * tr * radLen;
+ fvec EMASS(0.511 * 1e-3); // GeV
fvec BETA = P / sqrt(E2Corrected);
fvec GAMMA = sqrt(E2Corrected) / fMass;
// Calculate xi factor (KeV).
- fvec XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
+ 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 = 2. * EMASS * ETASQ;
- fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
- fvec EMAX = 1e6 * F1 / F2;
+ 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 * (1. - (BETA * BETA / 2.)) * 1e-12;
+ fvec DEDX2 = XI * EMAX * (fvec(1.) - (BETA * BETA / fvec(2.))) * fvec(1e-12);
- float P2 = P * P;
- fvec SDEDX = ((E2) *DEDX2) / (P2 * P2 * P2);
+ fvec P2 = P * P;
+ fvec SDEDX = (E2 * DEDX2) / (P2 * P2 * P2);
// T.C40 *= corr;
// T.C41 *= corr;
diff --git a/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.cxx b/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.cxx
index cfa04e1d88..26db240835 100644
--- a/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.cxx
+++ b/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.cxx
@@ -395,9 +395,9 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(const int NHits, nsL1vector<E
sHit.S2 = sHit.lx * sHit.ly;
sHit.S3 = sHit.lx * lr2;
sHit.S4 = sHit.ly * lr2;
- sHit.C = -lr * .5;
+ sHit.C = -lr * fvec(.5);
- const fvec w = iif(validHit, iif(lr > fvec(1.E-4), 1. / lr, fvec::One()), fvec::Zero());
+ const fvec w = iif(validHit, iif(lr > fvec(1.E-4), fvec(1.) / lr, fvec::One()), fvec::Zero());
const fvec w2 = w * w;
sHit.Cx = w * sHit.lx;
sHit.Cy = w * sHit.ly;
@@ -423,7 +423,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(const int NHits, nsL1vector<E
fvec tmp = S0 * S1 - S2 * S2;
// if( fabs(tmp) < 1.E-10 ) break; // CHECKME
- tmp = 0.5 / tmp;
+ tmp = fvec(0.5) / tmp;
X = (S3 * S1 - S4 * S2) * tmp;
Y = (S4 * S0 - S3 * S2) * tmp;
R2 = X * X + Y * Y;
@@ -431,7 +431,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(const int NHits, nsL1vector<E
const fvec Dcut = Dmax * Rejection; // cut for noise hits
// float RingCut = HitSize4 * R ; // closeness
- S0 = S1 = S2 = S3 = S4 = 0.0;
+ S0 = S1 = S2 = S3 = S4 = fvec::Zero();
Dmax = -1.;
for (THitIndex ih = 0; ih < MaxSearchAreaSize; ih++) {
@@ -446,7 +446,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(const int NHits, nsL1vector<E
const fvec dp = iif(sHit.on_ring, fvec(-1.), abs(sHit.C + sHit.Cx * X + sHit.Cy * Y));
Dmax = iif(((dp <= Dcut) & (dp > Dmax)), dp, Dmax);
- fvec w = iif((sHit.on_ring), 1. / (HitSizeSq_v + abs(dSq)), 1. / (1.e-5 + abs(dSq)));
+ fvec w = iif((sHit.on_ring), fvec(1.) / (HitSizeSq_v + abs(dSq)), fvec(1.) / (fvec(1.e-5) + abs(dSq)));
w = iif((dp <= Dcut) & validHit, w, fvec::Zero());
S0 += w * sHit.S0;
S1 += w * sHit.S1;
@@ -487,8 +487,8 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(const int NHits, nsL1vector<E
// if( fabs(s0) < 1.E-6 || fabs(s1) < 1.E-6 ) continue; // CHECKME
const fvec tmp = s1 * (S5 * S5 - NRingHits * S0) + s0 * s0;
const fvec A = ((S0 * S7 - S3 * S5) * s1 - s2 * s0) / tmp;
- Y = (s2 + s0 * A) / s1 * 0.5;
- X = (S3 + S5 * A - S2 * Y * 2) / S0 * 0.5;
+ Y = (s2 + s0 * A) / s1 * fvec(0.5);
+ X = (S3 + S5 * A - S2 * Y * 2) / S0 * fvec(0.5);
R2 = X * X + Y * Y - A;
// if( R2 < 0 ) continue; // will be rejected letter by R2 > R2Min
R = sqrt(R2);
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index 3d3c88e799..8ec8617c0a 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -48,10 +48,10 @@ using std::vector;
namespace NS_L1TrackFitter
{
- const fvec c_light = 0.000299792458, c_light_i = 1. / c_light;
- const fvec ZERO = 0.;
- const fvec ONE = 1.;
- const fvec vINF = 0.1;
+ const fvec c_light(0.000299792458), c_light_i = fvec(1.) / c_light;
+ const fvec ZERO = fvec(0.);
+ const fvec ONE = fvec(1.);
+ const fvec vINF = fvec(0.1);
} // namespace NS_L1TrackFitter
using namespace NS_L1TrackFitter;
@@ -498,7 +498,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<L1FieldRe
for (int iSt = nStations - 4; iSt >= 0; iSt--) {
- fvec w = iif(T.z > (zSta[iSt] + 2.5), fvec::One(), fvec::Zero());
+ fvec w = iif(T.z > (zSta[iSt] + fvec(2.5)), fvec::One(), fvec::Zero());
L1Extrapolate(T, zSta[iSt], T.qp, fld, &w);
if (iSt == NMvdStations - 1) {
@@ -518,14 +518,14 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<L1FieldRe
Double_t Cv[3] = {primVtx.GetCovMatrix()[0], primVtx.GetCovMatrix()[1], primVtx.GetCovMatrix()[2]};
- fvec dx = T.x - primVtx.GetRefX();
- fvec dy = T.y - primVtx.GetRefY();
+ fvec dx = T.x - fvec(primVtx.GetRefX());
+ fvec dy = T.y - fvec(primVtx.GetRefY());
fvec c[3] = {T.C00, T.C10, T.C11};
c[0] += fvec(Cv[0]);
c[1] += fvec(Cv[1]);
c[2] += fvec(Cv[2]);
fvec d = c[0] * c[2] - c[1] * c[1];
- fvec chi = sqrt(abs(0.5 * (dx * dx * c[0] - 2 * dx * dy * c[1] + dy * dy * c[2]) / d));
+ fvec chi = sqrt(abs(fvec(0.5) * (dx * dx * c[0] - fvec(2.) * dx * dy * c[1] + dy * dy * c[2]) / d));
chi.setZero(abs(d) < fvec(1.e-20));
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
diff --git a/reco/L1/vectors/L1vec.h b/reco/L1/vectors/L1vec.h
index 901ab2f7fa..d706949356 100644
--- a/reco/L1/vectors/L1vec.h
+++ b/reco/L1/vectors/L1vec.h
@@ -8,4 +8,6 @@
#include "vectors/L1vecVc.h"
//#include "vectors/L1vecPseudo.h"
+#include "std_alloc.h"
+
#endif
diff --git a/reco/L1/vectors/L1vecPseudo.h b/reco/L1/vectors/L1vecPseudo.h
index 9ea93f1604..ad0af0e0dc 100644
--- a/reco/L1/vectors/L1vecPseudo.h
+++ b/reco/L1/vectors/L1vecPseudo.h
@@ -119,7 +119,7 @@ public:
fscal v[fmask::Size];
- fvec() : fvec(0.f) {}
+ fvec() : fvec(0.) {}
fvec(const fvec& a)
{
@@ -235,22 +235,11 @@ public:
friend fvec operator-(const fvec& a) { return fvec(0) - a; }
friend fvec operator+(const fvec& a) { return a; }
- friend fvec operator+(const fvec& a, const fscal& b) { return a + fvec(b); }
- friend fvec operator-(const fvec& a, const fscal& b) { return a - fvec(b); }
- friend fvec operator*(const fvec& a, const fscal& b) { return a * fvec(b); }
- friend fvec operator/(const fvec& a, const fscal& b) { return a / fvec(b); }
- friend fvec operator+(const fscal& a, const fvec& b) { return fvec(a) + b; }
- friend fvec operator-(const fscal& a, const fvec& b) { return fvec(a) - b; }
- friend fvec operator*(const fscal& a, const fvec& b) { return fvec(a) * b; }
- friend fvec operator/(const fscal& a, const fvec& b) { return fvec(a) / b; }
+
friend void operator+=(fvec& a, const fvec& b) { a = a + b; }
friend void operator-=(fvec& a, const fvec& b) { a = a - b; }
friend void operator*=(fvec& a, const fvec& b) { a = a * b; }
friend void operator/=(fvec& a, const fvec& b) { a = a / b; }
- friend void operator+=(fvec& a, const fscal& b) { a = a + b; }
- friend void operator-=(fvec& a, const fscal& b) { a = a - b; }
- friend void operator*=(fvec& a, const fscal& b) { a = a * b; }
- friend void operator/=(fvec& a, const fscal& b) { a = a / b; }
friend std::ostream& operator<<(std::ostream& strm, const fvec& a)
{
@@ -272,9 +261,7 @@ public:
} __attribute__((aligned(16)));
-
#define _fvecalignment __attribute__((aligned(fvec::size() * sizeof(fscal))))
-#include "std_alloc.h"
#endif
diff --git a/reco/L1/vectors/L1vecVc.h b/reco/L1/vectors/L1vecVc.h
index b38e36bc45..639a221624 100644
--- a/reco/L1/vectors/L1vecVc.h
+++ b/reco/L1/vectors/L1vecVc.h
@@ -13,19 +13,4 @@ typedef Vc::float_m fmask;
#define _fvecalignment __attribute__((aligned(Vc::VectorAlignment)))
-inline fvec operator*(fscal a, const fvec& b) { return fvec(a) * b; }
-inline fvec operator*(const fvec& a, fscal b) { return a * fvec(b); }
-
-inline fvec operator/(fscal a, const fvec& b) { return fvec(a) / b; }
-inline fvec operator/(const fvec& a, fscal b) { return a / fvec(b); }
-
-inline fvec operator+(fscal a, const fvec& b) { return fvec(a) + b; }
-inline fvec operator+(const fvec& a, fscal b) { return a + fvec(b); }
-
-inline fvec operator-(fscal a, const fvec& b) { return fvec(a) - b; }
-inline fvec operator-(const fvec& a, fscal b) { return a - fvec(b); }
-
-#include "std_alloc.h"
-
-
#endif
--
GitLab