diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp
index 8f8774f07f0b476f2ede5aebd98a0d9e9058f048..af223589ab7e3e72fbe528a75ece76d2ecfd9356 100644
--- a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp
+++ b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp
@@ -54,13 +54,11 @@ CaloFutureSCorrection::CaloFutureSCorrection( const std::string& type, const std
 
 StatusCode CaloFutureSCorrection::finalize() {
   m_hypos.clear();
-  // finalize the base class
   return CaloFutureCorrectionBase::finalize();
 }
 // ============================================================================
 
 StatusCode CaloFutureSCorrection::initialize() {
-  // first initialize the base class
   StatusCode sc = CaloFutureCorrectionBase::initialize();
   if ( sc.isFailure() ) { return Error( "Unable initialize the base class CaloFutureCorrectionBase!", sc ); }
   if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) debug() << "Condition name : " << m_conditionName << endmsg;
@@ -78,11 +76,9 @@ StatusCode CaloFutureSCorrection::correct(
   }
   for ( auto* hypo : hypos ) {
 
-    // check the Hypo
     auto h = std::find( m_hypos.begin(), m_hypos.end(), hypo->hypothesis() );
     if ( m_hypos.end() == h ) return Error( "Invalid hypothesis!", StatusCode::SUCCESS );
 
-    // No correction for negative energy :
     if ( hypo->e() < 0. ) {
       ++m_counterSkipNegativeEnergyCorrection;
       continue;
@@ -113,13 +109,6 @@ StatusCode CaloFutureSCorrection::correct(
       continue;
     }
 
-    // int area = cellID.area(); // currently not used
-
-    // Somewhat inelegant way of data sharing between this const method process() and calcSCorrection(). Use of private
-    // fields of this CaloFutureSCorrection class for the same purpose would break the constness of the process()
-    // interface.
-    // Note, that calling private calcSCorrection(...,  &results) is thread-safe b/c that's a ptr to a stack variable.
-
     const LHCb::CaloPosition& position = MainCluster->position();
     const double              xBar     = position.x();
     const double              yBar     = position.y();
@@ -128,8 +117,6 @@ StatusCode CaloFutureSCorrection::correct(
     auto                      seedPos  = m_det->cellCenter( cellID );
     SCorrInputParams          params{cellID, seedPos, xBar, yBar, z};
 
-    // SCorrOutputParams results = {0, 0}; // output parameters: just d(Xhypo)/d(Xcluster) and d(Yhypo)/d(Ycluster)
-
     /** here all information is available
      *
      *  ( ) Ecal energy in 3x3     :   ( not used )
@@ -141,30 +128,28 @@ StatusCode CaloFutureSCorrection::correct(
      *  (7) Position of seed cell  :    seedPos
      */
 
-    // passing OUTPUT parameter &results as ptr to calcSCorrection() is thread-safe b/c that's a local stack variable
-    // calcSCorrection( xBar, yBar, xCor, yCor, params,
-    //&results ); // non-nullptr results means calculate and update results
-
     struct SCorrResults results = getSCorrections( params );
     double              xCor    = results.xCor;
     double              yCor    = results.yCor;
-    ; // corrected hypo position
-    const double& dXhy_dXcl = results.dXhy_dXcl;
-    const double& dYhy_dYcl = results.dYhy_dYcl;
+    // corrected hypo position
+    double dXhy_dXcl = results.dXhy_dXcl;
+    double dYhy_dYcl = results.dYhy_dYcl;
 
     // protection against unphysical d(Xhypo)/d(Xcluster) == 0 or d(Yhypo)/d(Ycluster) == 0
     if ( fabs( dXhy_dXcl ) < 1e-10 ) {
       warning() << "unphysical d(Xhypo)/d(Xcluster) = " << dXhy_dXcl << " reset to 1 as if Xhypo = Xcluster" << endmsg;
-      const_cast<double&>( dXhy_dXcl ) = 1.;
+      dXhy_dXcl = 1.;
     }
     if ( fabs( dYhy_dYcl ) < 1e-10 ) {
       warning() << "unphysical d(Yhypo)/d(Ycluster) = " << dYhy_dYcl << " reset to 1 as if Yhypo = Ycluster" << endmsg;
-      const_cast<double&>( dYhy_dYcl ) = 1.;
+      dYhy_dYcl = 1.;
     }
 
     // numeric partial derivatives w.r.t. X and Y, necessary to check after any change to the S-corrections
     if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) && m_correctCovariance ) {
-      debugDerivativesCalculation( xBar, yBar, xCor, yCor, dXhy_dXcl, dYhy_dYcl, params );
+      const double& dXhy_dXcl_test = dXhy_dXcl;
+      const double& dYhy_dYcl_test = dYhy_dYcl;
+      debugDerivativesCalculation( xBar, yBar, xCor, yCor, dXhy_dXcl_test, dYhy_dYcl_test, params );
     }
 
     const LHCb::CaloPosition* pos = hypo->position();
@@ -179,70 +164,59 @@ StatusCode CaloFutureSCorrection::correct(
       debug() << "xcel/ycel " << params.seedPos.x() << "/" << params.seedPos.y() << endmsg;
     }
 
-    // update position
-    LHCb::CaloPosition::Parameters& parameters = hypo->position()->parameters();
-    parameters( LHCb::CaloPosition::Index::X ) = xCor;
-    parameters( LHCb::CaloPosition::Index::Y ) = yCor;
+    updatePosition(xCor, yCor, hypo);
     m_counterDeltaX += xCor - xBar;
     m_counterDeltaY += yCor - yBar;
 
-    // update cov.m.: error propagation due to the S-correction
-    if ( m_correctCovariance ) {
-      LHCb::CaloPosition::Covariance& covariance = hypo->position()->covariance();
-
-      if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) { debug() << "before s-cor cov.m. = \n" << covariance << endmsg; }
-
-      // cov.m packing in double array[5] following ROOT::Math::SMatrix<double,3,3>::Array()
-      // for row/column indices (X:0, Y:1, E:2), see comments in CaloFutureECorrection::process()
-      double c1[6];
-
-      c1[0] = covariance( LHCb::CaloPosition::Index::X,
-                          LHCb::CaloPosition::Index::X ); // arr[0] not relying on LHCb::CaloPosition::Index::X == 0
-      c1[2] = covariance( LHCb::CaloPosition::Index::Y,
-                          LHCb::CaloPosition::Index::Y ); // arr[2] not relying on LHCb::CaloPosition::Index::Y == 1
-      c1[5] = covariance( LHCb::CaloPosition::Index::E,
-                          LHCb::CaloPosition::Index::E ); // arr[5] not relying on LHCb::CaloPosition::Index::E == 2
-      c1[1] = covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::Y ); // arr[1]
-      c1[3] = covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::E ); // arr[3]
-      c1[4] = covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::E ); // arr[4]
-
-      // cov1 = (J * cov0 * J^T) for the special case of diagonal Jacobian for (X,Y,E) -> (X1=X1(X), Y1=Y1(Y), E1=E)
-      c1[0] *= dXhy_dXcl * dXhy_dXcl;
-      c1[1] *= dXhy_dXcl * dYhy_dYcl;
-      c1[2] *= dYhy_dYcl * dYhy_dYcl;
-      c1[3] *= dXhy_dXcl;
-      c1[4] *= dYhy_dYcl;
-      // c1[5] remains unchanged (energy is not chaged by S-correction)
-
-      // alternatively, a code fragment for a general-form Jacobian (cf. a similar comment in
-      // CaloFutureECorrection::process()) TMatrixD jac(3, 3); // just a diagonal Jacobian in case of (X,Y,E) -> (X1(X),
-      // Y1(Y), E) transformation jac(0,0) = dXhy_dXcl; jac(1,1) = dYhy_dYcl; jac(2,2) = 1.; if ( msgLevel( MSG::DEBUG)
-      // ){ debug() << "s-cor jacobian = " << endmsg; jac.Print(); } TMarixDSym cov0(3) = ...          // to be
-      // initilized from hypo->position()->covariance() TMarixDSym cov1(3);               // resulting extrapolated
-      // cov.m. recalculate_cov(jac, cov0, cov1); // calculate:  cov1 = (J * cov0 * J^T)
-
-      // finally update CaloHypo::position()->covariance()
-      covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::X ) = c1[0]; // cov1(0,0);
-      covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::Y ) = c1[2]; // cov1(1,1);
-      covariance( LHCb::CaloPosition::Index::E, LHCb::CaloPosition::Index::E ) = c1[5]; // cov1(2,2);
-      covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::Y ) = c1[1]; // cov1(0,1);
-      covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::E ) = c1[3]; // cov1(0,2);
-      covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::E ) = c1[4]; // cov1(1,2);
-
-      // // DG: my little paranoia, should be always ok since Covariance is SMatrix<3,3,double> internally represented
-      // as double array[5] assert( covariance(LHCb::CaloFuturePosition::Index::X, LHCb::CaloFuturePosition::Index::Y)
-      // == covariance(LHCb::CaloFuturePosition::Index::Y, LHCb::CaloFuturePosition::Index::X)); assert(
-      // covariance(LHCb::CaloFuturePosition::Index::X, LHCb::CaloFuturePosition::Index::E) ==
-      // covariance(LHCb::CaloFuturePosition::Index::E, LHCb::CaloFuturePosition::Index::X)); assert(
-      // covariance(LHCb::CaloFuturePosition::Index::Y, LHCb::CaloFuturePosition::Index::E) ==
-      // covariance(LHCb::CaloFuturePosition::Index::E, LHCb::CaloFuturePosition::Index::Y));
-
-      if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) { debug() << "after s-cor cov.m. = \n" << covariance << endmsg; }
+    if ( m_correctCovariance) {
+      updateCovariance(dXhy_dXcl, dYhy_dYcl, hypo);
     }
   }
   return StatusCode::SUCCESS;
 }
 
+void CaloFutureSCorrection::updatePosition(double xCor, double yCor, LHCb::CaloHypo* hypo) const {
+  LHCb::CaloPosition::Parameters& parameters = hypo->position()->parameters();
+  parameters( LHCb::CaloPosition::Index::X ) = xCor;
+  parameters( LHCb::CaloPosition::Index::Y ) = yCor;
+}
+void CaloFutureSCorrection::updateCovariance(double dXhy_dXcl, double dYhy_dYcl, LHCb::CaloHypo* hypo) const {
+  LHCb::CaloPosition::Covariance& covariance = hypo->position()->covariance();
+
+  if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) { debug() << "before s-cor cov.m. = \n" << covariance << endmsg; }
+
+  // cov.m packing in double array[5] following ROOT::Math::SMatrix<double,3,3>::Array()
+  // for row/column indices (X:0, Y:1, E:2), see comments in CaloFutureECorrection::process()
+  double c1[6];
+
+  c1[0] = covariance( LHCb::CaloPosition::Index::X,
+                      LHCb::CaloPosition::Index::X ); // arr[0] not relying on LHCb::CaloPosition::Index::X == 0
+  c1[2] = covariance( LHCb::CaloPosition::Index::Y,
+                      LHCb::CaloPosition::Index::Y ); // arr[2] not relying on LHCb::CaloPosition::Index::Y == 1
+  c1[5] = covariance( LHCb::CaloPosition::Index::E,
+                      LHCb::CaloPosition::Index::E ); // arr[5] not relying on LHCb::CaloPosition::Index::E == 2
+  c1[1] = covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::Y ); // arr[1]
+  c1[3] = covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::E ); // arr[3]
+  c1[4] = covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::E ); // arr[4]
+
+  // cov1 = (J * cov0 * J^T) for the special case of diagonal Jacobian for (X,Y,E) -> (X1=X1(X), Y1=Y1(Y), E1=E)
+  c1[0] *= dXhy_dXcl * dXhy_dXcl;
+  c1[1] *= dXhy_dXcl * dYhy_dYcl;
+  c1[2] *= dYhy_dYcl * dYhy_dYcl;
+  c1[3] *= dXhy_dXcl;
+  c1[4] *= dYhy_dYcl;
+  // c1[5] remains unchanged (energy is not changed by S-correction)
+
+  covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::X ) = c1[0]; // cov1(0,0);
+  covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::Y ) = c1[2]; // cov1(1,1);
+  covariance( LHCb::CaloPosition::Index::E, LHCb::CaloPosition::Index::E ) = c1[5]; // cov1(2,2);
+  covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::Y ) = c1[1]; // cov1(0,1);
+  covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::E ) = c1[3]; // cov1(0,2);
+  covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::E ) = c1[4]; // cov1(1,2);
+
+  if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) { debug() << "after s-cor cov.m. = \n" << covariance << endmsg; }
+}
+
 // ============================================================================
 
 struct CaloFutureSCorrection::SCorrResults
@@ -348,37 +322,6 @@ CaloFutureSCorrection::getSCorrections( const struct SCorrInputParams& params )
    */
 
   if ( m_correctCovariance ) {
-    // // ---- calculation of numeric derivatives of individual correction functions, important for debugging in case of
-    // code changes --- debug() << "---------- numeric derivatives of individual S-correction functions ---------------"
-    // << endmsg; double tmpd = ( fabs(Asx0) > 1.e-5 ) ? Asx0*2.e-2 : 2.e-7; double dn_shapeX    = (
-    // getCorrection(CaloFutureCorrection::shapeX, cellID, Asx0 + tmpd, Asx0 + tmpd) - Asx1 )/tmpd; tmpd = ( fabs(Asy0)
-    // > 1.e-5 ) ? Asy0*2.e-2 : 2.e-7; double dn_shapeY    = ( getCorrection(CaloFutureCorrection::shapeY, cellID, Asy0
-    // + tmpd, Asy0 + tmpd) - Asy1 )/tmpd;
-    //
-    // double dn_angularX  = ( getCorrection(CaloFutureCorrection::angularX, cellID, thx*1.002, 0.) - daX )/thx/2e-3;
-    // double dn_angularY  = ( getCorrection(CaloFutureCorrection::angularY, cellID, thy*1.002, 0.) - daY )/thy/2e-3;
-    //
-    // tmpd = ( fabs(Asx2) > 1.e-5 ) ? Asx2*2.e-3 : 2.e-8;
-    // double dn_residualX = ( getCorrection((residualX_flag ? CaloFutureCorrection::residualX :
-    // CaloFutureCorrection::residual),
-    //                                                                                          cellID, Asx2 + tmpd, 0.)
-    //                                                                                          - dcX )/tmpd;
-    // tmpd = ( fabs(Asy2) > 1.e-5 ) ? Asy2*2.e-3 : 2.e-8;
-    // double dn_residualY = ( getCorrection((residualY_flag ? CaloFutureCorrection::residualY :
-    // CaloFutureCorrection::residual),
-    //                                                                                         cellID, Asy2 + tmpd, 0.)
-    //                                                                                         - dcY )/tmpd;
-    // tmpd = ( fabs(Asx3) > 1.e-5 ) ? Asx3*2.e-3 : 2.e-8;
-    // double dn_asymX     = (xBar < 0 ) ?
-    //   ( getCorrection(CaloFutureCorrection::asymM , cellID , Asx2 + tmpd , 0.) - ddcX )/tmpd :
-    //   ( getCorrection(CaloFutureCorrection::asymP , cellID , Asx2 + tmpd , 0.) - ddcX )/tmpd  ;
-    //
-    // tmpd = ( fabs(Asy3) > 1.e-5 ) ? Asy3*2.e-3 : 2.e-8;
-    // double dn_asymY     = (yBar < 0 ) ?
-    //   ( getCorrection(CaloFutureCorrection::asymM , cellID , Asy2 + tmpd , 0.) - ddcY )/tmpd :
-    //   ( getCorrection(CaloFutureCorrection::asymP , cellID , Asy2 + tmpd , 0.) - ddcY )/tmpd  ;
-    // //
-    // -------------------------------------------------------------------------------------------------------------------------------
 
     // calculation of the analytic derivatives:
     // NB: printouts comparing analytic calculations with numeric derivatives which are commented-out below
@@ -610,3 +553,50 @@ void CaloFutureSCorrection::debugDerivativesCalculation( const double& xBar, con
 
   return;
 }
+
+// WK: This was part of the comment - maybe it should be moved as kind of unit tests?
+// // ---- calculation of numeric derivatives of individual correction functions, important for debugging in case of
+// code changes --- debug() << "---------- numeric derivatives of individual S-correction functions ---------------"
+// << endmsg; double tmpd = ( fabs(Asx0) > 1.e-5 ) ? Asx0*2.e-2 : 2.e-7; double dn_shapeX    = (
+// getCorrection(CaloFutureCorrection::shapeX, cellID, Asx0 + tmpd, Asx0 + tmpd) - Asx1 )/tmpd; tmpd = ( fabs(Asy0)
+// > 1.e-5 ) ? Asy0*2.e-2 : 2.e-7; double dn_shapeY    = ( getCorrection(CaloFutureCorrection::shapeY, cellID, Asy0
+// + tmpd, Asy0 + tmpd) - Asy1 )/tmpd;
+//
+// double dn_angularX  = ( getCorrection(CaloFutureCorrection::angularX, cellID, thx*1.002, 0.) - daX )/thx/2e-3;
+// double dn_angularY  = ( getCorrection(CaloFutureCorrection::angularY, cellID, thy*1.002, 0.) - daY )/thy/2e-3;
+//
+// tmpd = ( fabs(Asx2) > 1.e-5 ) ? Asx2*2.e-3 : 2.e-8;
+// double dn_residualX = ( getCorrection((residualX_flag ? CaloFutureCorrection::residualX :
+// CaloFutureCorrection::residual),
+//                                                                                          cellID, Asx2 + tmpd, 0.)
+//                                                                                          - dcX )/tmpd;
+// tmpd = ( fabs(Asy2) > 1.e-5 ) ? Asy2*2.e-3 : 2.e-8;
+// double dn_residualY = ( getCorrection((residualY_flag ? CaloFutureCorrection::residualY :
+// CaloFutureCorrection::residual),
+//                                                                                         cellID, Asy2 + tmpd, 0.)
+//                                                                                         - dcY )/tmpd;
+// tmpd = ( fabs(Asx3) > 1.e-5 ) ? Asx3*2.e-3 : 2.e-8;
+// double dn_asymX     = (xBar < 0 ) ?
+//   ( getCorrection(CaloFutureCorrection::asymM , cellID , Asx2 + tmpd , 0.) - ddcX )/tmpd :
+//   ( getCorrection(CaloFutureCorrection::asymP , cellID , Asx2 + tmpd , 0.) - ddcX )/tmpd  ;
+//
+// tmpd = ( fabs(Asy3) > 1.e-5 ) ? Asy3*2.e-3 : 2.e-8;
+// double dn_asymY     = (yBar < 0 ) ?
+//   ( getCorrection(CaloFutureCorrection::asymM , cellID , Asy2 + tmpd , 0.) - ddcY )/tmpd :
+//   ( getCorrection(CaloFutureCorrection::asymP , cellID , Asy2 + tmpd , 0.) - ddcY )/tmpd  ;
+// //
+// -------------------------------------------------------------------------------------------------------------------------------
+//
+// // DG: my little paranoia, should be always ok since Covariance is SMatrix<3,3,double> internally represented
+// as double array[5] assert( covariance(LHCb::CaloFuturePosition::Index::X, LHCb::CaloFuturePosition::Index::Y)
+// == covariance(LHCb::CaloFuturePosition::Index::Y, LHCb::CaloFuturePosition::Index::X)); assert(
+// covariance(LHCb::CaloFuturePosition::Index::X, LHCb::CaloFuturePosition::Index::E) ==
+// covariance(LHCb::CaloFuturePosition::Index::E, LHCb::CaloFuturePosition::Index::X)); assert(
+// covariance(LHCb::CaloFuturePosition::Index::Y, LHCb::CaloFuturePosition::Index::E) ==
+// covariance(LHCb::CaloFuturePosition::Index::E, LHCb::CaloFuturePosition::Index::Y));
+// alternatively, a code fragment for a general-form Jacobian (cf. a similar comment in
+// CaloFutureECorrection::process()) TMatrixD jac(3, 3); // just a diagonal Jacobian in case of (X,Y,E) -> (X1(X),
+// Y1(Y), E) transformation jac(0,0) = dXhy_dXcl; jac(1,1) = dYhy_dYcl; jac(2,2) = 1.; if ( msgLevel( MSG::DEBUG)
+// ){ debug() << "s-cor jacobian = " << endmsg; jac.Print(); } TMarixDSym cov0(3) = ...          // to be
+// initilized from hypo->position()->covariance() TMarixDSym cov1(3);               // resulting extrapolated
+// cov.m. recalculate_cov(jac, cov0, cov1); // calculate:  cov1 = (J * cov0 * J^T)
diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h
index 85ee7f3bcbfe2e3a48d0ac6307db08d54c0c5b56..2616ecc44a5cb70c4e119d00d94a9f5e028a8b18 100644
--- a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h
+++ b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h
@@ -77,6 +77,9 @@ private:
 
   struct SCorrResults getSCorrections( const struct SCorrInputParams& params ) const;
 
+  void updateCovariance(double dXhy_dXcl, double dYhy_dYcl, LHCb::CaloHypo* hypo) const; 
+  void updatePosition(double xCor, double yCor, LHCb::CaloHypo* hypo) const; 
+
 private:
   using IncCounter    = Gaudi::Accumulators::Counter<>;
   using SCounter      = Gaudi::Accumulators::StatCounter<float>;