diff --git a/hist/hist/inc/TF1.h b/hist/hist/inc/TF1.h
index ab274bdb3b5ef..759457afe1948 100644
--- a/hist/hist/inc/TF1.h
+++ b/hist/hist/inc/TF1.h
@@ -29,9 +29,10 @@
 #include "TAttFill.h"
 #include "TAttMarker.h"
 #include "TMath.h"
-
+#include "Math/Math_vectypes.hxx"
 #include "Math/ParamFunctor.h"
 
+
 class TF1;
 class TH1;
 class TAxis;
@@ -250,6 +251,7 @@ class TF1 : public TNamed, public TAttLine, public TAttFill, public TAttMarker {
    Bool_t      fNormalized;  //Normalization option (false by default)
    Double_t    fNormIntegral;//Integral of the function before being normalized
    TF1FunctorPointer  *fFunctor = nullptr; //! Functor object to wrap any C++ callable object
+   TF1FunctionPointer *fFunctp = nullptr;  //! Pointer to vectorized function
    TFormula    *fFormula;    //Pointer to TFormula in case when user define formula
    TF1Parameters *fParams;   //Pointer to Function parameters object (exusts only for not-formula functions)
 
@@ -257,8 +259,8 @@ class TF1 : public TNamed, public TAttLine, public TAttFill, public TAttMarker {
 
    template<class T>
    struct TF1FunctionPointerImpl: TF1FunctionPointer {
-      TF1FunctionPointerImpl(const std::function<T(const T *f, const Double_t *param)> &&func): fimpl(func) {};
-      std::function<T(const T *f, const Double_t *param)> fimpl;
+      TF1FunctionPointerImpl(const std::function<T(const T *f, const Double_t *param)> &&func) : fImpl(func){};
+      std::function<T(const T *f, const Double_t *param)> fImpl;
    };
 
 
@@ -272,7 +274,6 @@ class TF1 : public TNamed, public TAttLine, public TAttFill, public TAttMarker {
 
 
 
-   TF1FunctionPointer *fFunctp = nullptr; //!Pointer to vectorized function
 
    static std::atomic<Bool_t> fgAbsValue;  //use absolute value of function when computing integral
    static Bool_t fgRejectPoint;  //True if point must be rejected in a fit
@@ -466,6 +467,9 @@ class TF1 : public TNamed, public TAttLine, public TAttFill, public TAttMarker {
    virtual Double_t EvalPar(const Double_t *x, const Double_t *params = 0);
    template<class T> T EvalPar(const T *x, const Double_t *params = 0);
    template<class T> T EvalParVec(const T *data, const Double_t *params = 0);
+#ifdef R__HAS_VECCORE
+   inline double EvalParVec(const Double_t *data, const Double_t *params);
+#endif
    virtual Double_t operator()(Double_t x, Double_t y = 0, Double_t z = 0, Double_t t = 0) const;
    virtual Double_t operator()(const Double_t *x, const Double_t *params = 0);
    template<class T> T operator()(const T *data, const Double_t *params);
@@ -814,19 +818,38 @@ T TF1::EvalPar(const T *x, const Double_t *params)
 template<class T>
 inline T TF1::EvalParVec(const T *data, const Double_t *params)
 {
-   // if (fType != 3) {
-   //    //This should throw an error
-   //    return TF1::EvalPar((double *) data, params);
-   // }
-   assert(fType == 3); 
-   if (fFunctor != nullptr)
+   assert(fType == 3);
+   if (!params) params = (Double_t *)fParams->GetParameters();
+   if (fFunctor)
       return ((TF1FunctorPointerImpl<T> *)fFunctor)->fImpl(data, params);
 
+   if (fFunctp)
+      return ((TF1FunctionPointerImpl<T> *)fFunctp)->fImpl(data, params);
+
    // this should throw an error
    // we nned to implement a vectorized GetSave(x)
-   return TMath::SignalingNaN(); 
+   return TMath::SignalingNaN();
+}
 
+#ifdef R__HAS_VECCORE
+inline double TF1::EvalParVec(const Double_t *data, const Double_t *params)
+{
+   assert(fType == 3);
+   ROOT::Double_v d, res;
+
+   d = *data;
+
+   if (fFunctor) {
+      res = ((TF1FunctorPointerImpl<ROOT::Double_v> *)fFunctor)->fImpl(&d, params);
+   } else if (fFunctp) {
+      res = ((TF1FunctionPointerImpl<ROOT::Double_v> *)fFunctp)->fImpl(&d, params);
+   } else {
+      //    res = GetSave(x);
+      return TMath::SignalingNaN();
+   }
+   return vecCore::Get<ROOT::Double_v>(res, 0);
 }
+#endif
 
 inline void TF1::SetRange(Double_t xmin, Double_t,  Double_t xmax, Double_t)
 {
diff --git a/hist/hist/src/HFitImpl.cxx b/hist/hist/src/HFitImpl.cxx
index 63883d1cabc40..e8b24c483eb2e 100644
--- a/hist/hist/src/HFitImpl.cxx
+++ b/hist/hist/src/HFitImpl.cxx
@@ -362,7 +362,7 @@ TFitResultPtr HFit::Fit(FitObject * h1, TF1 *f1 , Foption_t & fitOption , const
       fitConfig.SetWeightCorrection(weight);
       bool extended = ((fitOption.Like & 4 ) != 4 );
       //if (!extended) Info("HFitImpl","Do a not -extended binned fit");
-      fitok = fitter->LikelihoodFit(*fitdata, extended);
+      fitok = fitter->LikelihoodFit(*fitdata, extended, fitOption.ExecPolicy);
    }
    else{ // standard least square fit
       fitok = fitter->Fit(*fitdata, fitOption.ExecPolicy);
@@ -893,7 +893,7 @@ TFitResultPtr ROOT::Fit::UnBinFit(ROOT::Fit::UnBinData * data, TF1 * fitfunc, Fo
    bool extended = (fitOption.Like & 1) == 1;
 
    bool fitok = false;
-   fitok = fitter->LikelihoodFit(fitdata, extended);
+   fitok = fitter->LikelihoodFit(fitdata, extended, fitOption.ExecPolicy);
    if ( !fitok  && !fitOption.Quiet )
       Warning("UnBinFit","Abnormal termination of minimization.");
 
diff --git a/hist/hist/src/TF1.cxx b/hist/hist/src/TF1.cxx
index c09426a879d50..8d1f1b907f7a9 100644
--- a/hist/hist/src/TF1.cxx
+++ b/hist/hist/src/TF1.cxx
@@ -1259,7 +1259,7 @@ Double_t TF1::EvalPar(const Double_t *x, const Double_t *params)
    }
 
    if (fType == 3) {
-      if (fFunctor) {
+      if (fFunctor || fFunctp) {
          if (params) result =  EvalParVec(x, params);
          else result =  EvalParVec(x, (Double_t *) fParams->GetParameters());
       }
diff --git a/math/mathcore/inc/Fit/Chi2FCN.h b/math/mathcore/inc/Fit/Chi2FCN.h
index 3633bc096b7e9..80945b225f599 100644
--- a/math/mathcore/inc/Fit/Chi2FCN.h
+++ b/math/mathcore/inc/Fit/Chi2FCN.h
@@ -120,21 +120,13 @@ class Chi2FCN : public BasicFCN<DerivFunType, ModelFunType, BinData> {
    /// i-th chi-square residual
    virtual double DataElement(const double *x, unsigned int i, double *g) const {
       if (i==0) this->UpdateNCalls();
-#ifdef R__HAS_VECCORE
       return FitUtil::Evaluate<T>::EvalChi2Residual(BaseFCN::ModelFunction(), BaseFCN::Data(), x, i, g);
-#else
-      return FitUtil::EvaluateChi2Residual(BaseFCN::ModelFunction(), BaseFCN::Data(), x, i, g);
-#endif
    }
 
    // need to be virtual to be instantiated
    virtual void Gradient(const double *x, double *g) const {
       // evaluate the chi2 gradient
-#ifdef R__HAS_VECCORE
       FitUtil::Evaluate<T>::EvalChi2Gradient(BaseFCN::ModelFunction(), BaseFCN::Data(), x, g, fNEffPoints);
-#else
-      FitUtil::EvaluateChi2Gradient(BaseFCN::ModelFunction(), BaseFCN::Data(), x, g, fNEffPoints);
-#endif
    }
 
    /// get type of fit method function
@@ -154,15 +146,9 @@ class Chi2FCN : public BasicFCN<DerivFunType, ModelFunType, BinData> {
    virtual double DoEval (const double * x) const {
       this->UpdateNCalls();
       if (BaseFCN::Data().HaveCoordErrors() || BaseFCN::Data().HaveAsymErrors())
-#ifdef R__HAS_VECCORE
          return FitUtil::Evaluate<T>::EvalChi2Effective(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints);
       else
          return FitUtil::Evaluate<T>::EvalChi2(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints, fExecutionPolicy);
-#else
-        return  FitUtil::EvaluateChi2Effective(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints);
-      else
-        return FitUtil::EvaluateChi2(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints, fExecutionPolicy);
-#endif
    }
 
    // for derivatives
diff --git a/math/mathcore/inc/Fit/FitUtil.h b/math/mathcore/inc/Fit/FitUtil.h
index 2c0cbb8d148ee..277fd6e8ec65e 100644
--- a/math/mathcore/inc/Fit/FitUtil.h
+++ b/math/mathcore/inc/Fit/FitUtil.h
@@ -32,15 +32,17 @@
 #include "TError.h"
 #include "TSystem.h"
 
+// using parameter cache is not thread safe but needed for normalizing the functions
+#define USE_PARAMCACHE
+
 #ifdef R__HAS_VECCORE
 namespace vecCore{
-   //Auxiliar function. To be included in VecCore's new release
-   template <typename T> vecCore::Scalar<T> Reduce(const T &v)
-   {
-      vecCore::Scalar<T> sum{};
-      for (size_t i = 0; i < VectorSize<T>(); ++i)
-      sum += vecCore::Get(v, i);
-      return sum;
+template <class T>
+vecCore::Mask_v<T> Int2Mask(unsigned i)
+{
+   T x;
+   for (unsigned j = 0; j < vecCore::VectorSize<T>(); j++) vecCore::Set<T>(x, i, i + 1);
+   return vecCore::Mask_v<T>(x < i);
    }
 }
 #endif
@@ -335,9 +337,9 @@ namespace FitUtil {
 
    unsigned setAutomaticChunking(unsigned nEvents);
 
-#ifdef R__HAS_VECCORE
    template<class T>
-   struct Evaluate{
+   struct Evaluate {
+#ifdef R__HAS_VECCORE
       static double EvalChi2(const IModelFunctionTempl<T> &func, const BinData & data, const double * p, unsigned int &nPoints, const unsigned int &executionPolicy, unsigned nChunks = 0)
       {
          // evaluate the chi2 given a  vectorized function reference  , the data and returns the value and also in nPoints
@@ -422,6 +424,7 @@ namespace FitUtil {
             for (unsigned int i = 0; i < (data.Size() / vecSize); i++) {
                res += mapFunction(i);
             }
+
 #ifdef R__USE_IMT
          } else if (executionPolicy == ROOT::Fit::kMultithread) {
             auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(data.Size() / vecSize);
@@ -436,11 +439,11 @@ namespace FitUtil {
          }
          nPoints = n;
 
-#ifdef DEBUG
-         std::cout << "chi2 = " << chi2 << " n = " << nPoints  /*<< " rejected = " << nRejected */ << std::endl;
-#endif
+         // Last SIMD vector of elements (if padding needed)
+         vecCore::MaskedAssign(res, vecCore::Int2Mask<T>(data.Size() % vecSize),
+                               res + mapFunction(data.Size() / vecSize));
 
-         return vecCore::Reduce(res);
+         return vecCore::ReduceAdd(res);
       }
 
       static double EvalLogL(const IModelFunctionTempl<T> &func, const UnBinData & data, const double * const p, int iWeight,
@@ -478,7 +481,7 @@ namespace FitUtil {
                T xmin_v, xmax_v;
                vecCore::Load<T>(xmin_v, xmin.data());
                vecCore::Load<T>(xmax_v, xmax.data());
-               if (vecCore::Reduce(func(&xmin_v, p)) != 0 || vecCore::Reduce(func(&xmax_v, p)) != 0) {
+               if (vecCore::ReduceAdd(func(&xmin_v, p)) != 0 || vecCore::ReduceAdd(func(&xmax_v, p)) != 0) {
                   MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood", "A range has not been set and the function is not zero at +/- inf");
                   return 0;
                }
@@ -498,7 +501,7 @@ namespace FitUtil {
             if(data.NDim() > 1) {
                std::vector<T> x(data.NDim());
                for (unsigned int j = 0; j < data.NDim(); ++j)
-                  vecCore::Load<T>(x[j], data.GetCoordComponent(i, j));
+                  vecCore::Load<T>(x[j], data.GetCoordComponent(i * vecSize, j));
 #ifdef USE_PARAMCACHE
                fval = func(x.data());
 #else
@@ -507,7 +510,7 @@ namespace FitUtil {
                // one -dim case
             } else {
                T x;
-               vecCore::Load<T>(x, data.GetCoordComponent(i, 0));
+               vecCore::Load<T>(x, data.GetCoordComponent(i * vecSize, 0));
 #ifdef USE_PARAMCACHE
                fval = func(&x);
 #else
@@ -612,7 +615,7 @@ namespace FitUtil {
                   T xmin_v, xmax_v;
                   vecCore::Load<T>(xmin_v, xmin.data());
                   vecCore::Load<T>(xmax_v, xmax.data());
-                  if (vecCore::Reduce(func(&xmin_v, p)) != 0 || vecCore::Reduce(func(&xmax_v, p)) != 0) {
+                  if (vecCore::ReduceAdd(func(&xmin_v, p)) != 0 || vecCore::ReduceAdd(func(&xmax_v, p)) != 0) {
                      MATH_ERROR_MSG("FitUtil::EvaluateLogLikelihood", "A range has not been set and the function is not zero at +/- inf");
                      return 0;
                   }
@@ -647,7 +650,7 @@ namespace FitUtil {
       }
 
       static double EvalPoissonLogL(const IModelFunctionTempl<T> &func, const BinData &data, const double *p, int iWeight, bool extended,
-                                    unsigned int &nPoints, const unsigned int &executionPolicy, unsigned nChunks = 0)
+                                    unsigned int, const unsigned int &executionPolicy, unsigned nChunks = 0)
       {
          // evaluate the Poisson Log Likelihood
          // for binned likelihood fits
@@ -664,10 +667,8 @@ namespace FitUtil {
          // iWeight = 2 ==> logL = Sum( w*w * f(x_i) )
          //
 
-         unsigned int n = data.Size();
-
 #ifdef USE_PARAMCACHE
-         (const_cast<IModelFunction &>(func)).SetParameters(p);
+         (const_cast<IModelFunctionTempl<T> &>(func)).SetParameters(p);
 #endif
          auto vecSize = vecCore::VectorSize<T>();
          // get fit option and check case of using integral of bins
@@ -685,7 +686,7 @@ namespace FitUtil {
             if (data.NDim() > 1) {
                std::vector<T> x(data.NDim());
                for (unsigned int j = 0; j < data.NDim(); ++j)
-                  vecCore::Load<T>(x[j], data.GetCoordComponent(i, j));
+                  vecCore::Load<T>(x[j], data.GetCoordComponent(i * vecSize, j));
 #ifdef USE_PARAMCACHE
                fval = func(x.data());
 #else
@@ -694,7 +695,7 @@ namespace FitUtil {
                // one -dim case
             } else {
                T x;
-               vecCore::Load<T>(x, data.GetCoordComponent(i, 0));
+               vecCore::Load<T>(x, data.GetCoordComponent(i * vecSize, 0));
 #ifdef USE_PARAMCACHE
                fval = func(&x);
 #else
@@ -704,8 +705,7 @@ namespace FitUtil {
 
             // EvalLog protects against 0 values of fval but don't want to add in the -log sum
             // negative values of fval
-            auto m = vecCore::Mask_v<T>(fval < 0.0);
-            vecCore::MaskedAssign<T>(fval, m, 0.0);
+            vecCore::MaskedAssign<T>(fval, fval < 0.0, 0.0);
 
             T nloglike{}; // negative loglikelihood
 
@@ -715,27 +715,21 @@ namespace FitUtil {
                // can apply correction only when y is not zero otherwise weight is undefined
                // (in case of weighted likelihood I don't care about the constant term due to
                // the saturated model)
-               if (y != 0) {
-                  T error{};
-                  auto pError = data.ErrorPtr(i * vecSize);
-                  std::vector<double> ones{1, 1, 1, 1};
-                  pError = (pError != nullptr) ? pError : &ones.front();
-                  vecCore::Load<T>(error, pError);
-
-                  T weight = (error * error) / y; // this is the bin effective weight
+               auto m = vecCore::Mask_v<T>(y != 0.0);
+               if (!vecCore::MaskFull(m)) {
+                  T error = 1;
+                  if (data.GetErrorType() != ROOT::Fit::BinData::ErrorType::kNoError)
+                     vecCore::Load<T>(error, data.ErrorPtr(i * vecSize));
+                  T weight;
+                  vecCore::MaskedAssign<T>(weight, y != 0, (error * error) / y);
                   if (extended) {
                      nloglike = fval * weight;
                      // wTot  += weight;
                      // w2Tot += weight*weight;
                   }
-                  nloglike -= weight * y * ROOT::Math::Util::EvalLog(fval);
+                  vecCore::MaskedAssign<T>(nloglike, y != 0, nloglike - weight * y * ROOT::Math::Util::EvalLog(fval));
                }
 
-               //  need to compute total weight and weight-square
-               // if (extended ) {
-               //    nuTot += fval;
-               // }
-
             } else {
                // standard case no weights or iWeight=1
                // this is needed for Poisson likelihood (which are extened and not for multinomial)
@@ -776,7 +770,11 @@ namespace FitUtil {
                "Execution policy unknown. Avalaible choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
          }
 
-         return vecCore::Reduce(res);
+         // Last padded SIMD vector of elements
+         vecCore::MaskedAssign(res, vecCore::Int2Mask<T>(data.Size() % vecSize),
+                               res + mapFunction(data.Size() / vecSize));
+
+         return vecCore::ReduceAdd(res);
       }
 
       static double EvalChi2Effective(const IModelFunctionTempl<T> &, const BinData &, const double *, unsigned int &)
@@ -795,10 +793,23 @@ namespace FitUtil {
          Error("FitUtil::Evaluate<T>::EvalChi2Residual", "The vectorized evaluation of the Chi2 with the ith residual is still not supported");
          return -1.;
       }
+
+      /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
+      /// and its gradient
+static double EvalPoissonBinPdf(const IModelFunctionTempl<T> &, const BinData &, const double *, unsigned int , double * ) {
+         Error("FitUtil::Evaluate<T>::EvaluatePoissonBinPdf", "The vectorized evaluation of the BinnedLikelihood fit evaluated point by point is still not supported");
+         return -1.;
+      }
+
+static void EvalPoissonLogLGradient(const IModelFunctionTempl<T> &, const BinData &, const double *, double * ) {
+         Error("FitUtil::Evaluate<T>::EvaluatePoissonLogLGradient", "The vectorized evaluation of the BinnedLikelihood fit evaluated point by point is still not supported");
+      }
    };
 
    template<>
    struct Evaluate<double>{
+#endif
+
       static double EvalChi2(const IModelFunction & func, const BinData & data, const double * p, unsigned int &nPoints, const unsigned int &executionPolicy, unsigned nChunks = 0)
       {
          // evaluate the chi2 given a  function reference, the data and returns the value and also in nPoints
@@ -807,7 +818,7 @@ namespace FitUtil {
          // optionally the integral of function in the bin is used
          return FitUtil::EvaluateChi2(func, data, p, nPoints, executionPolicy, nChunks);
       }
-      
+
       static double EvalLogL(const IModelFunctionTempl<double> &func, const UnBinData & data, const double * p, int iWeight,
       bool extended, unsigned int &nPoints, const unsigned int &executionPolicy, unsigned nChunks = 0)
       {
@@ -832,8 +843,17 @@ namespace FitUtil {
       {
          return FitUtil::EvaluateChi2Residual(func, data, p, i, g);
       }
+
+      /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
+      /// and its gradient
+      static double EvalPoissonBinPdf(const IModelFunctionTempl<double> &func, const BinData & data, const double *p, unsigned int i, double *g ) {
+         return FitUtil::EvaluatePoissonBinPdf(func, data, p, i, g);
+      }
+
+static void EvalPoissonLogLGradient(const IModelFunctionTempl<double> &func, const BinData &data, const double *p, double *g) {
+         FitUtil::EvaluatePoissonLogLGradient(func, data, p, g);
+      }
    };
-#endif
 
 } // end namespace FitUtil
 
diff --git a/math/mathcore/inc/Fit/Fitter.h b/math/mathcore/inc/Fit/Fitter.h
index fdf3209d2ba83..87da7861844ce 100644
--- a/math/mathcore/inc/Fit/Fitter.h
+++ b/math/mathcore/inc/Fit/Fitter.h
@@ -79,8 +79,12 @@ class Fitter {
 public:
 
    typedef ROOT::Math::IParamMultiFunction                 IModelFunction;
+   template <class T>
+   using IModelFunctionTempl =                             ROOT::Math::IParamMultiFunctionTempl<T>;
 #ifdef R__HAS_VECCORE
    typedef ROOT::Math::IParametricFunctionMultiDimTempl<ROOT::Double_v>  IModelFunction_v;
+#else
+   typedef ROOT::Math::IParamMultiFunction                 IModelFunction_v;
 #endif
    typedef ROOT::Math::IParamMultiGradFunction             IGradModelFunction;
    typedef ROOT::Math::IParamFunction                      IModel1DFunction;
@@ -169,7 +173,7 @@ class Fitter {
       SetData(data);
       return DoBinnedLikelihoodFit(extended, executionPolicy);
    }
-   
+
    bool LikelihoodFit(const std::shared_ptr<BinData> &data, bool extended = true,
                       ROOT::Fit::ExecutionPolicy executionPolicy = ROOT::Fit::kSerial) {
       SetData(data);
@@ -452,9 +456,10 @@ class Fitter {
       fData = std::shared_ptr<FitData>(const_cast<FitData*>(&data),DummyDeleter<FitData>());
    }
    // set data and function without cloning them
-   void SetFunctionAndData(const IModelFunction & func, const FitData & data) {
+   template <class T>
+   void SetFunctionAndData(const IModelFunctionTempl<T> & func, const FitData & data) {
       SetData(data);
-      fFunc = std::shared_ptr<IModelFunction>(const_cast<IModelFunction*>(&func),DummyDeleter<IModelFunction>());
+      fFunc = std::shared_ptr<IModelFunctionTempl<T>>(const_cast<IModelFunctionTempl<T>*>(&func),DummyDeleter<IModelFunctionTempl<T>>());
    }
 
    //set data for the fit using a shared ptr
@@ -487,11 +492,9 @@ class Fitter {
    int fDataSize;  // size of data sets (need for Fumili or LM fitters)
 
    FitConfig fConfig;       // fitter configuration (options and parameter settings)
-#ifdef R__HAS_VECCORE
+
    std::shared_ptr<IModelFunction_v> fFunc_v;  //! copy of the fitted  function containing on output the fit result
-#else
-    std::shared_ptr<IModelFunction> fFunc_v;   //dummy for when VecCore not available. Keeps the code cleaner.
-#endif
+
    std::shared_ptr<IModelFunction> fFunc;  //! copy of the fitted  function containing on output the fit result
 
    std::shared_ptr<ROOT::Fit::FitResult>  fResult;  //! pointer to the object containing the result of the fit
diff --git a/math/mathcore/inc/Fit/LogLikelihoodFCN.h b/math/mathcore/inc/Fit/LogLikelihoodFCN.h
index e14015e7b984f..47308a6b8b24e 100644
--- a/math/mathcore/inc/Fit/LogLikelihoodFCN.h
+++ b/math/mathcore/inc/Fit/LogLikelihoodFCN.h
@@ -154,12 +154,7 @@ class LogLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,UnBinData>  {
     */
    virtual double DoEval (const double * x) const {
       this->UpdateNCalls();
-
-#ifdef R__HAS_VECCORE
       return FitUtil::Evaluate<T>::EvalLogL(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fWeight, fIsExtended, fNEffPoints, fExecutionPolicy);
-#else
-      return FitUtil::EvaluateLogL(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fWeight, fIsExtended, fNEffPoints, fExecutionPolicy);
-#endif
    }
 
    // for derivatives
diff --git a/math/mathcore/inc/Fit/PoissonLikelihoodFCN.h b/math/mathcore/inc/Fit/PoissonLikelihoodFCN.h
index e67665d9bb207..99c78ab1dcb15 100644
--- a/math/mathcore/inc/Fit/PoissonLikelihoodFCN.h
+++ b/math/mathcore/inc/Fit/PoissonLikelihoodFCN.h
@@ -54,8 +54,8 @@ class PoissonLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,BinData>
    typedef  ::ROOT::Math::BasicFitMethodFunction<DerivFunType> BaseObjFunction;
    typedef typename  BaseObjFunction::BaseFunction BaseFunction;
 
-   typedef  ::ROOT::Math::IParamMultiFunction IModelFunction;
-
+   typedef ::ROOT::Math::IParamMultiFunctionTempl<T> IModelFunction;
+   typedef typename BaseObjFunction::Type_t Type_t;
 
    /**
       Constructor from unbin data set and model function (pdf)
@@ -122,13 +122,13 @@ class PoissonLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,BinData>
    /// i-th likelihood element and its gradient
    virtual double DataElement(const double * x, unsigned int i, double * g) const {
       if (i==0) this->UpdateNCalls();
-      return FitUtil::EvaluatePoissonBinPdf(BaseFCN::ModelFunction(), BaseFCN::Data(), x, i, g);
+      return FitUtil::Evaluate<typename BaseFCN::T>::EvalPoissonBinPdf(BaseFCN::ModelFunction(), BaseFCN::Data(), x, i, g);
    }
 
    /// evaluate gradient
    virtual void Gradient(const double *x, double *g) const {
       // evaluate the chi2 gradient
-      FitUtil::EvaluatePoissonLogLGradient(BaseFCN::ModelFunction(), BaseFCN::Data(), x, g );
+      FitUtil::Evaluate<typename BaseFCN::T>::EvalPoissonLogLGradient(BaseFCN::ModelFunction(), BaseFCN::Data(), x, g);
    }
 
    /// get type of fit method function
@@ -161,13 +161,8 @@ class PoissonLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,BinData>
     */
    virtual double DoEval (const double * x) const {
       this->UpdateNCalls();
-#ifdef R__HAS_VECCORE
       return FitUtil::Evaluate<T>::EvalPoissonLogL(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fWeight, fIsExtended,
                                                    fNEffPoints, fExecutionPolicy);
-#else
-      return FitUtil::EvaluatePoissonLogL(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fWeight, fIsExtended,
-                                          fNEffPoints, fExecutionPolicy);
-#endif
    }
 
    // for derivatives
diff --git a/math/mathcore/src/BinData.cxx b/math/mathcore/src/BinData.cxx
index fccebc37a79a3..9df9d76f01b46 100644
--- a/math/mathcore/src/BinData.cxx
+++ b/math/mathcore/src/BinData.cxx
@@ -627,8 +627,14 @@ namespace ROOT {
 
     void BinData::InitDataVector ()
     {
-      fData.resize( fMaxPoints );
-      fDataPtr = &fData.front();
+#ifdef R__HAS_VECCORE
+       // Add padding to be a multiple of SIMD vector size and help looping
+       auto extraP = vecCore::VectorSize<ROOT::Double_v>() - ((fMaxPoints) % vecCore::VectorSize<ROOT::Double_v>());
+       fData.resize(fMaxPoints + extraP);
+#else
+       fData.resize(fMaxPoints);
+#endif
+       fDataPtr = &fData.front();
     }
 
     void BinData::InitializeErrors()
diff --git a/math/mathcore/src/FitUtil.cxx b/math/mathcore/src/FitUtil.cxx
index 42bbd32d363a4..07c4cc01286c6 100644
--- a/math/mathcore/src/FitUtil.cxx
+++ b/math/mathcore/src/FitUtil.cxx
@@ -39,9 +39,6 @@
 #include <iostream>
 #endif
 
-// using parameter cache is not thread safe but needed for normalizing the functions
-#define USE_PARAMCACHE
-
 //  need to implement integral option
 
 namespace ROOT {
diff --git a/math/mathcore/src/Fitter.cxx b/math/mathcore/src/Fitter.cxx
index 1d2bd47611f33..6c3625f97ce62 100644
--- a/math/mathcore/src/Fitter.cxx
+++ b/math/mathcore/src/Fitter.cxx
@@ -339,12 +339,10 @@ bool Fitter::DoLeastSquareFit(ROOT::Fit::ExecutionPolicy executionPolicy) {
       if (!fFunc_v ) {
          MATH_ERROR_MSG("Fitter::DoLeastSquareFit","model function is not set");
          return false;
-#ifdef R__HAS_VECCORE
       } else{
          Chi2FCN<BaseFunc, IModelFunction_v> chi2(data, fFunc_v, executionPolicy);
          fFitType = chi2.Type();
          return DoMinimization (chi2);
-#endif
      }
    } else {
 
@@ -387,10 +385,11 @@ bool Fitter::DoBinnedLikelihoodFit(bool extended, ROOT::Fit::ExecutionPolicy exe
    bool  useWeight = fConfig.UseWeightCorrection();
 
    // check function
-   if (!fFunc) {
-      MATH_ERROR_MSG("Fitter::DoBinnedLikelihoodFit","model function is not set");
-      return false;
-   }
+   if (!fFunc)
+      if (!fFunc_v) {
+         MATH_ERROR_MSG("Fitter::DoBinnedLikelihoodFit", "model function is not set");
+         return false;
+      }
 
    // logl fit (error should be 0.5) set if different than default values (of 1)
    if (fConfig.MinimizerOptions().ErrorDef() == gDefaultErrorDef ) {
@@ -406,21 +405,35 @@ bool Fitter::DoBinnedLikelihoodFit(bool extended, ROOT::Fit::ExecutionPolicy exe
    fBinFit = true;
    fDataSize = data->Size();
 
-   // create a chi2 function to be used for the equivalent chi-square
-   Chi2FCN<BaseFunc> chi2(data,fFunc);
 
    if (!fUseGradient) {
       // do minimization without using the gradient
-      PoissonLikelihoodFCN<BaseFunc> logl(data, fFunc, useWeight, extended, executionPolicy);
-      fFitType = logl.Type();
-      // do minimization
-      if (!DoMinimization (logl, &chi2) ) return false;
-      if (useWeight) {
-         logl.UseSumOfWeightSquare();
-         if (!ApplyWeightCorrection(logl) ) return false;
+      if (fFunc_v) {
+         // create a chi2 function to be used for the equivalent chi-square
+         Chi2FCN<BaseFunc, IModelFunction_v> chi2(data, fFunc_v);
+         PoissonLikelihoodFCN<BaseFunc, IModelFunction_v> logl(data, fFunc_v, useWeight, extended, executionPolicy);
+         fFitType = logl.Type();
+         // do minimization
+         if (!DoMinimization(logl, &chi2)) return false;
+         if (useWeight) {
+            logl.UseSumOfWeightSquare();
+            if (!ApplyWeightCorrection(logl)) return false;
+         }
+      } else {
+         // create a chi2 function to be used for the equivalent chi-square
+         Chi2FCN<BaseFunc> chi2(data, fFunc);
+         PoissonLikelihoodFCN<BaseFunc> logl(data, fFunc, useWeight, extended, executionPolicy);
+         fFitType = logl.Type();
+         // do minimization
+         if (!DoMinimization(logl, &chi2)) return false;
+         if (useWeight) {
+            logl.UseSumOfWeightSquare();
+            if (!ApplyWeightCorrection(logl)) return false;
+         }
       }
-   }
-   else {
+   } else {
+      // create a chi2 function to be used for the equivalent chi-square
+      Chi2FCN<BaseFunc> chi2(data, fFunc);
       if (fConfig.MinimizerOptions().PrintLevel() > 0)
          MATH_INFO_MSG("Fitter::DoLikelihoodFit","use gradient from model function");
       // check if fFunc provides gradient
@@ -482,7 +495,6 @@ bool Fitter::DoUnbinnedLikelihoodFit(bool extended, ROOT::Fit::ExecutionPolicy e
 
    if (!fUseGradient) {
       // do minimization without using the gradient
-#ifdef R__HAS_VECCORE
      if (fFunc_v ){
        LogLikelihoodFCN<BaseFunc, IModelFunction_v> logl(data, fFunc_v, useWeight, extended, executionPolicy);
        fFitType = logl.Type();
@@ -492,10 +504,8 @@ bool Fitter::DoUnbinnedLikelihoodFit(bool extended, ROOT::Fit::ExecutionPolicy e
           if (!ApplyWeightCorrection(logl) ) return false;
         }
         return true;
-     }
-     else{         
-#endif
-       LogLikelihoodFCN<BaseFunc> logl(data, fFunc, useWeight, extended, executionPolicy);
+     } else {
+        LogLikelihoodFCN<BaseFunc> logl(data, fFunc, useWeight, extended, executionPolicy);
 
         fFitType = logl.Type();
         if (!DoMinimization (logl) ) return false;
@@ -504,9 +514,7 @@ bool Fitter::DoUnbinnedLikelihoodFit(bool extended, ROOT::Fit::ExecutionPolicy e
           if (!ApplyWeightCorrection(logl) ) return false;
         }
         return true;
-#ifdef R__HAS_VECCORE
      }
-#endif
    } else {
       // use gradient : check if fFunc provides gradient
       if (fConfig.MinimizerOptions().PrintLevel() > 0)
diff --git a/math/mathcore/test/fit/testBinnedFitExecPolicy.cxx b/math/mathcore/test/fit/testBinnedFitExecPolicy.cxx
index daa068d0671a8..ee78645322af4 100644
--- a/math/mathcore/test/fit/testBinnedFitExecPolicy.cxx
+++ b/math/mathcore/test/fit/testBinnedFitExecPolicy.cxx
@@ -5,6 +5,7 @@
 #include "Math/IParamFunction.h"
 #include "TFitResult.h"
 #include "TError.h"
+#include "Math/MinimizerOptions.h"
 
 int compareResult(double v1, double v2, std::string s = "", double tol = 0.01)
 {
@@ -26,74 +27,90 @@ int main()
 {
    TF1 *f = new TF1("fvCore", func<double>, 100, 200, 4);
    f->SetParameters(1, 1000, 7.5, 1.5);
-   TH1D h1f("h1f", "Test random numbers", 128000, 100, 200);
+
+   // NBins not multiple of SIMD vector size, testing padding
+   TH1D h1f("h1f", "Test random numbers", 12801, 100, 200);
    gRandom->SetSeed(1);
    h1f.FillRandom("fvCore", 1000000);
 
+   std::cout << "\n **FIT: Chi2 **\n\n";
    auto r1 = h1f.Fit(f, "S");
    if ((Int_t)r1 != 0) {
-      Error("testChi2ExecPolicy", "Sequential Chi2 Fit failed!");
+      Error("testBinnedFitExecPolicy", "Sequential Chi2 Fit failed!");
       return -1;
    }
 
+   std::cout << "\n **FIT: Binned Likelihood **\n\n";
    auto rL1 = h1f.Fit(f, "S L");
-   if ((Int_t)r1 != 0) {
-      Error("testChi2ExecPolicy", "Sequential  Binned Likelihood Fit failed!");
+   if ((Int_t)rL1 != 0) {
+      Error("testBinnedFitExecPolicy", "Sequential  Binned Likelihood Fit failed!");
       return -1;
    }
 
 #ifdef R__USE_IMT
+   std::cout << "\n **FIT: Multithreaded Chi2 **\n\n";
+   f->SetParameters(1, 1000, 7.5, 1.5);
    auto r2 = h1f.Fit(f, "MULTITHREAD S");
    if ((Int_t)r2 != 0) {
-      Error("testChi2ExecPolicy", "Multithreaded Chi2 Fit failed!");
+      Error("testBinnedFitExecPolicy", "Multithreaded Chi2 Fit failed!");
       return -1;
    } else {
-      compareResult(r2->Chi2(), r1->Chi2(), "Mutithreaded Chi2 Fit: ");
+      compareResult(r2->MinFcnValue(), r1->MinFcnValue(), "Mutithreaded Chi2 Fit: ");
    }
 
+   std::cout << "\n **FIT: Multithreaded Binned Likelihood **\n\n";
+   f->SetParameters(1, 1000, 7.5, 1.5);
    auto rL2 = h1f.Fit(f, "MULTITHREAD S L");
    if ((Int_t)rL2 != 0) {
-      Error("testChi2ExecPolicy", "Multithreaded Binned Likelihood Fit failed!");
+      Error("testBinnedFitExecPolicy", "Multithreaded Binned Likelihood Fit failed!");
       return -1;
    } else {
-      compareResult(rL2->Chi2(), rL1->Chi2(), "Mutithreaded Binned Likelihood Fit (PoissonLogL): ");
+      compareResult(rL2->MinFcnValue(), rL1->MinFcnValue(), "Mutithreaded Binned Likelihood Fit (PoissonLogL): ");
    }
 #endif
 
 #ifdef R__HAS_VECCORE
    TF1 *fvecCore = new TF1("fvCore", func<ROOT::Double_v>, 100, 200, 4);
+
+   std::cout << "\n **FIT: Vectorized Chi2 **\n\n";
    fvecCore->SetParameters(1, 1000, 7.5, 1.5);
    auto r3 = h1f.Fit(fvecCore, "S");
    if ((Int_t)r3 != 0) {
-      Error("testChi2ExecPolicy", "Vectorized Chi2 Fit failed!");
+      Error("testBinnedFitExecPolicy", "Vectorized Chi2 Fit failed!");
       return -1;
    } else {
-      compareResult(r3->Chi2(), r1->Chi2(), "Vectorized Chi2 Fit: ");
+      compareResult(r3->MinFcnValue(), r1->MinFcnValue(), "Vectorized Chi2 Fit: ");
    }
 
+   std::cout << "\n **FIT: Vectorized Binned Likelihood **\n\n";
+   fvecCore->SetParameters(1, 1000, 7.5, 1.5);
    auto rL3 = h1f.Fit(fvecCore, "S L");
    if ((Int_t)rL3 != 0) {
-      Error("testChi2ExecPolicy", "Vectorized Binned Likelihood Fit failed!");
+      Error("testBinnedFitExecPolicy", "Vectorized Binned Likelihood Fit failed!");
       return -1;
    } else {
-      compareResult(rL3->Chi2(), rL1->Chi2(), "Vectorized Binned Likelihood Fit (PoissonLogL) Fit: ");
+      compareResult(rL3->MinFcnValue(), rL1->MinFcnValue(), "Vectorized Binned Likelihood Fit (PoissonLogL) Fit: ");
    }
 
 #ifdef R__USE_IMT
+   std::cout << "\n **FIT: Mutithreaded vectorized Chi2 **\n\n";
    auto r4 = h1f.Fit(fvecCore, "MULTITHREAD S");
    if ((Int_t)r4 != 0) {
-      Error("testChi2ExecPolicy", "Multithreaded vectorized Fit failed!");
+      Error("testBinnedFitExecPolicy", "Mutithreaded vectorized Chi2 Fit failed!");
       return -1;
    } else {
-      compareResult(r4->Chi2(), r1->Chi2(), "Mutithreaded vectorized Chi2 Fit: ");
+      compareResult(r4->MinFcnValue(), r1->MinFcnValue(), "Mutithreaded vectorized Chi2 Fit: ");
    }
 
+   std::cout << "\n **FIT: Multithreaded and vectorized Binned Likelihood **\n\n";
+   fvecCore->SetParameters(1, 1000, 7.5, 1.5);
    auto rL4 = h1f.Fit(fvecCore, "MULTITHREAD S L");
    if ((Int_t)rL4 != 0) {
-      Error("testChi2ExecPolicy", "Multithreaded Binned Likelihood vectorized Fit failed!");
+      Error("testBinnedFitExecPolicy", "Multithreaded Binned Likelihood vectorized Fit failed!");
       return -1;
    } else {
-      compareResult(rL4->Chi2(), rL1->Chi2(), "Mutithreaded vectorized Binned Likelihood Fit (PoissonLogL) Fit: ");
+      compareResult(rL4->MinFcnValue(), rL1->MinFcnValue(),
+                    "Mutithreaded vectorized Binned Likelihood Fit (PoissonLogL) Fit: ");
    }
 
 #endif