diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 64466fa09f2111c2fe04a7b1f6e52ae780dbc1f5..e31c98fd5f4a75a79cdd5af06113eb2bbabe7069 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -164,7 +164,7 @@ inline void L1Algo::f11(  /// input 1st stage of singlet search
     fvec zl          = zPos_l[i1_V];
     fvec& time       = HitTime_l[i1_V];
     fvec& timeEr     = HitTimeEr[i1_V];
-    const fvec& dzli = 1. / (zl - targZ);
+    const fvec dzli  = 1. / (zl - targZ);
 
     fvec dx1, dy1, dxy1 = 0;
 
@@ -174,8 +174,8 @@ inline void L1Algo::f11(  /// input 1st stage of singlet search
 
     StripsToCoor(u, v, xl, yl, stal);
 
-    const fvec& tx = (xl - targX) * dzli;
-    const fvec& ty = (yl - targY) * dzli;
+    const fvec tx = (xl - targX) * dzli;
+    const fvec ty = (yl - targY) * dzli;
 
     // estimate field for singlet creation
     int istac = istal / 2;  // "center" station
@@ -402,12 +402,12 @@ inline void L1Algo::f20(  // input
 
     const int n2Saved = n2;
 
-    const fvec& Pick_m22 =
+    const fvec Pick_m22 =
       (DOUBLET_CHI2_CUT - T1.chi2);  // if make it bigger the found hits will be rejected later because of the chi2 cut.
     // Pick_m22 is not used, search for mean squared, 2nd version
 
     // -- collect possible doublets --
-    const fscal& iz        = 1 / T1.z[i1_4];
+    const fscal iz         = 1 / T1.z[i1_4];
     const float& timeError = T1.C55[i1_4];
     const float& time      = T1.t[i1_4];
 
@@ -439,7 +439,7 @@ inline void L1Algo::f20(  // input
       if ((Event[i1_V][i1_4] != hitm.n)) continue;
 #endif
       // - check whether hit belong to the window ( track position +\- errors ) -
-      const fscal& zm   = hitm.Z();
+      const fscal zm    = hitm.Z();
       L1TrackPar T1_new = T1;
       fvec dz           = fvec(zm) - T1.z;
 
@@ -467,7 +467,7 @@ inline void L1Algo::f20(  // input
 
       StripsToCoor(hitm.U(), hitm.V(), xm, ym, stam);
 
-      const fscal& dY = ym[i1_4] - y[i1_4];
+      const fscal dY = ym[i1_4] - y[i1_4];
 
       if (dY * dY > dy_est2 && dY < 0) continue;
 
@@ -483,7 +483,7 @@ inline void L1Algo::f20(  // input
         dx_est2 = Pick_m22[i1_4] * fabs(C00[i1_4] + dxm[0] * dxm[0]);
       }
 
-      const fscal& dX = xm[i1_4] - x[i1_4];
+      const fscal dX = xm[i1_4] - x[i1_4];
       if (dX * dX > dx_est2) continue;
 
       // check chi2
@@ -688,7 +688,7 @@ inline void L1Algo::f30(  // input
           continue;
         if (T2.C00[i2_4] < 0 || T2.C11[i2_4] < 0 || T2.C22[i2_4] < 0 || T2.C33[i2_4] < 0) continue;
 
-        const fvec& Pick_r22   = (TRIPLET_CHI2_CUT - T2.chi2);
+        const fvec Pick_r22    = (TRIPLET_CHI2_CUT - T2.chi2);
         const float& timeError = T2.C55[i2_4];
         const float& time      = T2.t[i2_4];
 
@@ -697,7 +697,7 @@ inline void L1Algo::f30(  // input
 #ifdef DO_NOT_SELECT_TRIPLETS
         if (isec == TRACKS_FROM_TRIPLETS_ITERATION) Pick_r22 = Pick_r2 + 1;
 #endif
-        const fscal& iz = 1 / T2.z[i2_4];
+        const fscal iz = 1 / T2.z[i2_4];
         L1HitAreaTime area(vGridTime[&star - vStations], T2.x[i2_4] * iz, T2.y[i2_4] * iz,
                            (sqrt(Pick_r22 * (T2.C00 + stam.XYInfo.C00)) + MaxDZ * fabs(T2.tx))[i2_4] * iz,
                            (sqrt(Pick_r22 * (T2.C11 + stam.XYInfo.C11)) + MaxDZ * fabs(T2.ty))[i2_4] * iz, time,
@@ -723,8 +723,8 @@ inline void L1Algo::f30(  // input
 #ifdef USE_EVENT_NUMBER
           if ((T2.n[i2_4] != hitr.n)) continue;
 #endif
-          const fscal& zr = hitr.Z();
-          //  const fscal& yr = hitr.Y();
+          const fscal zr = hitr.Z();
+          //  const fscal yr = hitr.Y();
 
           fvec xr, yr = 0;
 
@@ -756,8 +756,8 @@ inline void L1Algo::f30(  // input
             dy_est2 = (Pick_r22[i2_4] * (fabs(C11[i2_4] + dyr[0] * dyr[0])));
           }
 
-          const fscal& dY  = yr[i2_4] - y[i2_4];
-          const fscal& dY2 = dY * dY;
+          const fscal dY  = yr[i2_4] - y[i2_4];
+          const fscal dY2 = dY * dY;
           if (dY2 > dy_est2 && dY2 < 0) continue;  // if (yr < y_minus_new[i2_4]) continue;
 
           // check upper boundary
@@ -774,7 +774,7 @@ inline void L1Algo::f30(  // input
             dx_est2 = (Pick_r22[i2_4] * (fabs(C00[i2_4] + dxr[0] * dxr[0])));
           }
 
-          const fscal& dX = xr[i2_4] - x[i2_4];
+          const fscal dX = xr[i2_4] - x[i2_4];
           if (dX * dX > dx_est2) continue;
           // check chi2  // not effective
           fvec C10;
@@ -2636,8 +2636,8 @@ inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best
       if ((new_trip.GetMHit() != curr_trip->GetRHit())) continue;
       if ((new_trip.GetLHit() != curr_trip->GetMHit())) continue;
 
-      const fscal& qp1 = curr_trip->GetQp();
-      const fscal& qp2 = new_trip.GetQp();
+      const fscal qp1  = curr_trip->GetQp();
+      const fscal qp2  = new_trip.GetQp();
       fscal dqp        = fabs(qp1 - qp2);
       fscal Cqp        = curr_trip->Cqp;
       Cqp += new_trip.Cqp;
diff --git a/reco/L1/L1Algo/L1FitMaterial.h b/reco/L1/L1Algo/L1FitMaterial.h
index 2483c44802b741ea619214a6b31bd735e9863c69..6c0a2d193d075266048a787980d8c6eaa88cde0e 100644
--- a/reco/L1/L1Algo/L1FitMaterial.h
+++ b/reco/L1/L1Algo/L1FitMaterial.h
@@ -128,20 +128,20 @@ inline fvec L1Fit::ApproximateBetheBloch(const fvec& bg2)
   // 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;
+  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;
 
   const float mK   = 0.307075e-3f;  // [GeV*cm^2/g]
   const 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 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
+  const fvec maxT  = _2me * bg2;  // neglecting the electron mass
 
   //*** Density effect
   fvec d2(0.f);
@@ -151,7 +151,7 @@ inline fvec L1Fit::ApproximateBetheBloch(const fvec& bg2)
   fvec init = x > x1;
 
   d2            = fvec(init & (lhwI + x - 0.5f));
-  const fvec& r = (x1 - x) / (x1 - x0);
+  const fvec r  = (x1 - x) / (x1 - x0);
   init          = (x > x0) & (x1 > x);
   d2            = fvec(init & (lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r)) + fvec((!init) & d2);
 
@@ -185,11 +185,11 @@ inline fvec L1Fit::ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const
   const float mK   = 0.307075e-3f;  // [GeV*cm^2/g]
   const 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 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
+  const fvec maxT  = _2me * bg2;  // neglecting the electron mass
 
   //*** Density effect
   fvec d2(0.f);
@@ -199,7 +199,7 @@ inline fvec L1Fit::ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const
   fvec init = x > x1;
 
   d2            = fvec(init & (lhwI + x - 0.5f));
-  const fvec& r = (x1 - x) / (x1 - x0);
+  const fvec r  = (x1 - x) / (x1 - x0);
   init          = (x > x0) & (x1 > x);
   d2            = fvec(init & (lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r)) + fvec((!init) & d2);
 
@@ -224,12 +224,12 @@ inline float L1Fit::CalcQpAfterEloss(float qp, float eloss, float mass2)
 
 inline 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;
+  const fvec p2 = 1.f / (T.qp * T.qp);
+  const fvec E2 = fMass2 + p2;
 
   // fvec qp =T.qp;
 
-  const fvec& bethe = ApproximateBetheBloch(p2 / fMass2);
+  const fvec bethe = ApproximateBetheBloch(p2 / fMass2);
 
   // fvec bethe2 = BetheBlochAl(qp[0]);
 
@@ -242,7 +242,7 @@ inline void L1Fit::EnergyLossCorrection(L1TrackPar& T, const fvec& radThick, fve
   // dE = fabs(dE_);//dE2;//bethe * radThick*tr * 9.34961f*2.33f ;
 
 
-  const fvec& E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+  const fvec E2Corrected  = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
   fvec corr               = sqrt(p2 / (E2Corrected - fMass2));
   fvec init               = fvec(!(corr == corr)) | fvec(w < 1);
   corr                    = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
@@ -262,8 +262,8 @@ inline void L1Fit::EnergyLossCorrection(L1TrackPar& T, const fvec& radThick, fve
 
 inline void L1Fit::EnergyLossCorrectionAl(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;
+  const fvec p2 = 1.f / (T.qp * T.qp);
+  const fvec E2 = fMass2 + p2;
 
   fvec qp = T.qp;
 
@@ -283,9 +283,9 @@ inline void L1Fit::EnergyLossCorrectionAl(L1TrackPar& T, const fvec& radThick, f
   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, rho, 0.20, 3.00, i, atomicZ / atomicA);
 
-  // const fvec& bethe = ApproximateBetheBloch( p2/fMass2 );
+  // const fvec bethe = ApproximateBetheBloch( p2/fMass2 );
 
   // fvec bethe2 = BetheBlochAl(qp[0]);
 
@@ -296,7 +296,7 @@ inline void L1Fit::EnergyLossCorrectionAl(L1TrackPar& T, const fvec& radThick, f
   fvec dE = bethe * radThick * tr * radLen * rho;
 
 
-  const fvec& E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+  const fvec E2Corrected  = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
   fvec corr               = sqrt(p2 / (E2Corrected - fMass2));
   fvec init               = fvec(!(corr == corr)) | fvec(w < 1);
   corr                    = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
@@ -352,8 +352,8 @@ inline void L1Fit::EnergyLossCorrectionAl(L1TrackPar& T, const fvec& radThick, f
 
 inline void L1Fit::EnergyLossCorrectionCarbon(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;
+  const fvec p2 = 1.f / (T.qp * T.qp);
+  const fvec E2 = fMass2 + p2;
 
   fvec qp = T.qp;
 
@@ -368,10 +368,10 @@ inline void L1Fit::EnergyLossCorrectionCarbon(L1TrackPar& T, const fvec& radThic
   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, rho, 0.20, 3.00, i, atomicZ / atomicA);
 
 
-  // const fvec& bethe = ApproximateBetheBloch( p2/fMass2 );
+  // const fvec bethe = ApproximateBetheBloch( p2/fMass2 );
 
   //  fvec bethe2 = BetheBlochCarbon(qp[0]);
 
@@ -382,7 +382,7 @@ inline void L1Fit::EnergyLossCorrectionCarbon(L1TrackPar& T, const fvec& radThic
   fvec dE = bethe * radThick * tr * rho * radLen;
   //   fvec dE2 = bethe2 * radThick*tr*2.70f* 18.76f;
 
-  const fvec& E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+  const fvec E2Corrected  = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
   fvec corr               = sqrt(p2 / (E2Corrected - fMass2));
   fvec init               = fvec(!(corr == corr)) | fvec(w < 1);
   corr                    = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
@@ -440,8 +440,8 @@ inline void L1Fit::EnergyLossCorrectionCarbon(L1TrackPar& T, const fvec& radThic
 
 inline void L1Fit::EnergyLossCorrectionIron(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;
+  const fvec p2 = 1.f / (T.qp * T.qp);
+  const fvec E2 = fMass2 + p2;
 
   fvec qp = T.qp;
 
@@ -455,7 +455,7 @@ inline void L1Fit::EnergyLossCorrectionIron(L1TrackPar& T, const fvec& radThick,
   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, rho, 0.20, 3.00, i, atomicZ / atomicA);
 
   // fvec bethe2 = BetheBlochIron(T.qp[0]);
 
@@ -466,7 +466,7 @@ inline void L1Fit::EnergyLossCorrectionIron(L1TrackPar& T, const fvec& radThick,
   fvec dE = bethe * radThick * tr * radLen * rho;
   //fvec dE2 = bethe2 * radThick*tr * 1.758f*2.33;
 
-  const fvec& E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+  const fvec E2Corrected  = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
   fvec corr               = sqrt(p2 / (E2Corrected - fMass2));
   fvec init               = fvec(!(corr == corr)) | fvec(w < 1);
   corr                    = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx
index eb00bb70cc0dc829b9c0b16b22a446e53ef72f5e..c20d8138169124e7d0293a7e1596b0c2d21035dd 100644
--- a/reco/L1/L1Algo/L1TrackParFit.cxx
+++ b/reco/L1/L1Algo/L1TrackParFit.cxx
@@ -901,17 +901,16 @@ void L1TrackParFit::L1AddMaterial(L1MaterialInfo& info, fvec qp0, fvec w)
 
 void L1TrackParFit::EnergyLossCorrection(const fvec& radThick, fvec& qp0, fvec direction, fvec w)
 {
-  const fvec& p2 = 1.f / (fqp * fqp);
-  const fvec& E2 = fMass2 + p2;
+  const fvec p2 = 1.f / (fqp * fqp);
+  const fvec E2 = fMass2 + p2;
 
-  const fvec& bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2);
+  const fvec bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2);
 
   fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
 
-  const fvec& dE = bethe * radThick * tr * 2.33f * 9.34961f;
+  const fvec dE = bethe * radThick * tr * 2.33f * 9.34961f;
 
-  const fvec& E2Corrected =
-    (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+  const fvec E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
   fvec corr = sqrt(p2 / (E2Corrected - fMass2));
   fvec init = fvec(!(corr == corr)) | fvec(w < 1);
   corr      = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
@@ -928,8 +927,8 @@ void L1TrackParFit::EnergyLossCorrection(const fvec& radThick, fvec& qp0, fvec d
 
 void L1TrackParFit::EnergyLossCorrectionIron(const fvec& radThick, fvec& qp0, fvec direction, fvec w)
 {
-  const fvec& p2 = 1.f / (qp0 * qp0);
-  const fvec& E2 = fMass2 + p2;
+  const fvec p2 = 1.f / (qp0 * qp0);
+  const fvec E2 = fMass2 + p2;
 
   int atomicZ   = 26;
   float atomicA = 55.845f;
@@ -942,15 +941,14 @@ void L1TrackParFit::EnergyLossCorrectionIron(const fvec& radThick, fvec& qp0, fv
   else
     i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
 
-  const fvec& bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+  const fvec bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
 
   fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
 
 
   fvec dE = bethe * radThick * tr * radLen * rho;
 
-  const fvec& E2Corrected =
-    (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+  const fvec E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
   fvec corr = sqrt(p2 / (E2Corrected - fMass2));
   fvec init = fvec(!(corr == corr)) | fvec(w < 1);
   corr      = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
@@ -995,8 +993,8 @@ void L1TrackParFit::EnergyLossCorrectionIron(const fvec& radThick, fvec& qp0, fv
 
 void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& radThick, fvec& qp0, fvec direction, fvec w)
 {
-  const fvec& p2 = 1.f / (qp0 * qp0);
-  const fvec& E2 = fMass2 + p2;
+  const fvec p2 = 1.f / (qp0 * qp0);
+  const fvec E2 = fMass2 + p2;
 
   int atomicZ   = 6;
   float atomicA = 12.011f;
@@ -1009,15 +1007,14 @@ void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& radThick, fvec& qp0,
   else
     i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
 
-  const fvec& bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+  const fvec bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
 
   fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
 
 
   fvec dE = bethe * radThick * tr * radLen * rho;
 
-  const fvec& E2Corrected =
-    (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+  const fvec E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
   fvec corr = sqrt(p2 / (E2Corrected - fMass2));
   fvec init = fvec(!(corr == corr)) | fvec(w < 1);
   corr      = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
@@ -1062,8 +1059,8 @@ void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& radThick, fvec& qp0,
 
 void L1TrackParFit::EnergyLossCorrectionAl(const fvec& radThick, fvec& qp0, fvec direction, fvec w)
 {
-  const fvec& p2 = 1.f / (qp0 * qp0);
-  const fvec& E2 = fMass2 + p2;
+  const fvec p2 = 1.f / (qp0 * qp0);
+  const fvec E2 = fMass2 + p2;
 
   int atomicZ   = 13;
   float atomicA = 26.981f;
@@ -1076,15 +1073,14 @@ void L1TrackParFit::EnergyLossCorrectionAl(const fvec& radThick, fvec& qp0, fvec
   else
     i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
 
-  const fvec& bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+  const fvec bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
 
   fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
 
 
   fvec dE = bethe * radThick * tr * radLen * rho;
 
-  const fvec& E2Corrected =
-    (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+  const fvec E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
   fvec corr = sqrt(p2 / (E2Corrected - fMass2));
   fvec init = fvec(!(corr == corr)) | fvec(w < 1);
   corr      = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));