diff --git a/roofit/roofit/inc/RooDstD0BG.h b/roofit/roofit/inc/RooDstD0BG.h
index c53e475d545bab..c9024facdb2204 100644
--- a/roofit/roofit/inc/RooDstD0BG.h
+++ b/roofit/roofit/inc/RooDstD0BG.h
@@ -45,6 +45,8 @@ class RooDstD0BG : public RooAbsPdf {
   RooRealProxy C,A,B ;
 
   Double_t evaluate() const;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
+
 
 private:
 
diff --git a/roofit/roofit/inc/RooLognormal.h b/roofit/roofit/inc/RooLognormal.h
index 114f3b6529cc2f..12a2d458a693f9 100644
--- a/roofit/roofit/inc/RooLognormal.h
+++ b/roofit/roofit/inc/RooLognormal.h
@@ -38,7 +38,8 @@ class RooLognormal : public RooAbsPdf {
   RooRealProxy k ;
 
   Double_t evaluate() const ;
-
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
+  
 private:
 
   ClassDef(RooLognormal,1) // log-normal PDF
diff --git a/roofit/roofit/inc/RooNovosibirsk.h b/roofit/roofit/inc/RooNovosibirsk.h
index 1a7259c1a26c61..6d19c0bbe15ea2 100644
--- a/roofit/roofit/inc/RooNovosibirsk.h
+++ b/roofit/roofit/inc/RooNovosibirsk.h
@@ -45,10 +45,12 @@ class RooNovosibirsk : public RooAbsPdf {
 
 protected:
   RooRealProxy x;
-  RooRealProxy width;
   RooRealProxy peak;
+  RooRealProxy width;
   RooRealProxy tail;
   Double_t evaluate() const;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
+
 
 private:
   ClassDef(RooNovosibirsk,1) // Novosibirsk PDF
diff --git a/roofit/roofit/inc/RooUniform.h b/roofit/roofit/inc/RooUniform.h
index 530d3af8b2567a..08f4f344c0b59d 100644
--- a/roofit/roofit/inc/RooUniform.h
+++ b/roofit/roofit/inc/RooUniform.h
@@ -40,6 +40,9 @@ class RooUniform : public RooAbsPdf {
   RooListProxy x ;
 
   Double_t evaluate() const ;
+  inline RooSpan<double> evaluateBatch(std::size_t, std::size_t ) const {
+    return {};
+  }
 
 private:
 
diff --git a/roofit/roofit/inc/RooVoigtian.h b/roofit/roofit/inc/RooVoigtian.h
index 8f0d259eaa1c9d..449dc39e6a126a 100644
--- a/roofit/roofit/inc/RooVoigtian.h
+++ b/roofit/roofit/inc/RooVoigtian.h
@@ -47,12 +47,13 @@ class RooVoigtian : public RooAbsPdf {
   RooRealProxy sigma ;
 
   Double_t evaluate() const ;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
+
 
 private:
 
-  Double_t _invRootPi;
   Bool_t _doFast;
-  ClassDef(RooVoigtian,1) // Voigtian PDF (Gauss (x) BreitWigner)
+  ClassDef(RooVoigtian,2) // Voigtian PDF (Gauss (x) BreitWigner)
 };
 
 #endif
diff --git a/roofit/roofit/src/RooDstD0BG.cxx b/roofit/roofit/src/RooDstD0BG.cxx
index f941c673183075..b0cf96cd013747 100644
--- a/roofit/roofit/src/RooDstD0BG.cxx
+++ b/roofit/roofit/src/RooDstD0BG.cxx
@@ -23,19 +23,18 @@ Special p.d.f shape that can be used to model the background of
 D*-D0 mass difference distributions
 **/
 
-#include "RooFit.h"
-
-#include "Riostream.h"
-#include "Riostream.h"
-#include <math.h>
-#include "TMath.h"
-
 #include "RooDstD0BG.h"
+#include "RooFit.h"
 #include "RooAbsReal.h"
 #include "RooRealVar.h"
 #include "RooIntegrator1D.h"
 #include "RooAbsFunc.h"
+#include "RooVDTHeaders.h"
+#include "BatchHelpers.h"
 
+#include "TMath.h"
+
+#include <cmath>
 using namespace std;
 
 ClassImp(RooDstD0BG);
@@ -74,6 +73,53 @@ Double_t RooDstD0BG::evaluate() const
   return (val > 0 ? val : 0) ;
 }
 
+namespace {
+//Author: Emmanouil Michalainas, CERN 9 SEPTEMBER 2019 
+
+template<class Tdm, class Tdm0, class TC, class TA, class TB>
+void compute(	size_t batchSize, double * __restrict output,
+              Tdm DM, Tdm0 DM0, TC C, TA A, TB B)
+{
+  for (size_t i=0; i<batchSize; i++) {
+    const double ratio = DM[i] / DM0[i];
+    const double arg1 = (DM0[i]-DM[i]) / C[i];
+    const double arg2 = A[i]*_rf_fast_log(ratio);
+    output[i] = (1 -_rf_fast_exp(arg1)) * _rf_fast_exp(arg2) +B[i]*(ratio-1);
+  }
+  
+  for (size_t i=0; i<batchSize; i++) {
+    if (output[i]<0) output[i] = 0;
+  }
+}
+};
+
+RooSpan<double> RooDstD0BG::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
+  using namespace BatchHelpers;
+
+  EvaluateInfo info = getInfo( {&dm, &dm0, &C, &A, &B}, begin, batchSize );
+  if (info.nBatches == 0) {
+    return {};
+  }
+  auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
+  auto dmData = dm.getValBatch(begin, info.size);
+
+  if (info.nBatches==1 && !dmData.empty()) {
+    compute(info.size, output.data(), dmData.data(),
+    BracketAdapter<double> (dm0),
+    BracketAdapter<double> (C),
+    BracketAdapter<double> (A),
+    BracketAdapter<double> (B));
+  }
+  else {
+    compute(info.size, output.data(),
+    BracketAdapterWithMask (dm,dm.getValBatch(begin,info.size)),
+    BracketAdapterWithMask (dm0,dm0.getValBatch(begin,info.size)),
+    BracketAdapterWithMask (C,C.getValBatch(begin,info.size)),
+    BracketAdapterWithMask (A,A.getValBatch(begin,info.size)),
+    BracketAdapterWithMask (B,B.getValBatch(begin,info.size)));
+  }
+  return output;
+}
 
 ////////////////////////////////////////////////////////////////////////////////
 /// if (matchArgs(allVars,analVars,dm)) return 1 ;
diff --git a/roofit/roofit/src/RooLognormal.cxx b/roofit/roofit/src/RooLognormal.cxx
index 495599c1c0b860..e8c8b5c7361932 100644
--- a/roofit/roofit/src/RooLognormal.cxx
+++ b/roofit/roofit/src/RooLognormal.cxx
@@ -23,25 +23,22 @@ The parameterization here is physics driven and differs from the ROOT::Math::log
   - x0 = 0
 **/
 
-#include "RooFit.h"
-
-#include "Riostream.h"
-#include "Riostream.h"
-#include <math.h>
-
 #include "RooLognormal.h"
+#include "RooFit.h"
 #include "RooAbsReal.h"
 #include "RooRealVar.h"
 #include "RooRandom.h"
 #include "RooMath.h"
-#include "TMath.h"
+#include "RooVDTHeaders.h"
+#include "RooHelpers.h"
+#include "BatchHelpers.h"
 
+#include "TMath.h"
 #include <Math/SpecFuncMathCore.h>
 #include <Math/PdfFuncMathCore.h>
 #include <Math/ProbFuncMathCore.h>
 
-#include "TError.h"
-
+#include <cmath>
 using namespace std;
 
 ClassImp(RooLognormal);
@@ -56,6 +53,14 @@ RooLognormal::RooLognormal(const char *name, const char *title,
   m0("m0","m0",this,_m0),
   k("k","k",this,_k)
 {
+    RooHelpers::checkRangeOfParameters(this, {&_x, &_m0, &_k}, 0.);
+    
+    auto par = dynamic_cast<const RooAbsRealLValue*>(&_k);
+    if (par && par->getMin()<=1 && par->getMax()>=1 ) {
+      oocoutE(this, InputArguments) << "The parameter '" << par->GetName() << "' with range [" << par->getMin("") << ", "
+          << par->getMax() << "] of the " << this->IsA()->GetName() << " '" << this->GetName()
+          << "' can reach the unsafe value 1.0 " << ". Advise to limit its range." << std::endl;
+    }
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -74,15 +79,75 @@ RooLognormal::RooLognormal(const RooLognormal& other, const char* name) :
 
 Double_t RooLognormal::evaluate() const
 {
-  Double_t xv = x;
   Double_t ln_k = TMath::Abs(TMath::Log(k));
   Double_t ln_m0 = TMath::Log(m0);
-  Double_t x0 = 0;
 
-  Double_t ret = ROOT::Math::lognormal_pdf(xv,ln_m0,ln_k,x0);
+  Double_t ret = ROOT::Math::lognormal_pdf(x,ln_m0,ln_k);
   return ret ;
 }
 
+////////////////////////////////////////////////////////////////////////////////
+
+namespace {
+//Author: Emmanouil Michalainas, CERN 10 September 2019
+
+template<class Tx, class Tm0, class Tk>
+void compute(	size_t batchSize,
+              double * __restrict output,
+              Tx X, Tm0 M0, Tk K)
+{
+  const double rootOf2pi = std::sqrt(2 * M_PI);
+  for (size_t i=0; i<batchSize; i++) {
+    double lnxOverM0 = _rf_fast_log(X[i]/M0[i]);
+    double lnk = _rf_fast_log(K[i]);
+    if (lnk<0) lnk = -lnk;
+    double arg = lnxOverM0/lnk;
+    arg *= -0.5*arg;
+    output[i] = _rf_fast_exp(arg) / (X[i]*lnk*rootOf2pi);
+  }
+}
+};
+
+RooSpan<double> RooLognormal::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
+  using namespace BatchHelpers;
+  auto xData = x.getValBatch(begin, batchSize);
+  auto m0Data = m0.getValBatch(begin, batchSize);
+  auto kData = k.getValBatch(begin, batchSize);
+  const bool batchX = !xData.empty();
+  const bool batchM0 = !m0Data.empty();
+  const bool batchK = !kData.empty();
+
+  if (!batchX && !batchM0 && !batchK) {
+    return {};
+  }
+  batchSize = findSize({ xData, m0Data, kData });
+  auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
+
+  if (batchX && !batchM0 && !batchK ) {
+    compute(batchSize, output.data(), xData, BracketAdapter<double>(m0), BracketAdapter<double>(k));
+  }
+  else if (!batchX && batchM0 && !batchK ) {
+    compute(batchSize, output.data(), BracketAdapter<double>(x), m0Data, BracketAdapter<double>(k));
+  }
+  else if (batchX && batchM0 && !batchK ) {
+    compute(batchSize, output.data(), xData, m0Data, BracketAdapter<double>(k));
+  }
+  else if (!batchX && !batchM0 && batchK ) {
+    compute(batchSize, output.data(), BracketAdapter<double>(x), BracketAdapter<double>(m0), kData);
+  }
+  else if (batchX && !batchM0 && batchK ) {
+    compute(batchSize, output.data(), xData, BracketAdapter<double>(m0), kData);
+  }
+  else if (!batchX && batchM0 && batchK ) {
+    compute(batchSize, output.data(), BracketAdapter<double>(x), m0Data, kData);
+  }
+  else if (batchX && batchM0 && batchK ) {
+    compute(batchSize, output.data(), xData, m0Data, kData);
+  }
+  return output;
+}
+
+
 ////////////////////////////////////////////////////////////////////////////////
 
 Int_t RooLognormal::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
diff --git a/roofit/roofit/src/RooNovosibirsk.cxx b/roofit/roofit/src/RooNovosibirsk.cxx
index 17c25483f07452..13c0b8cd4ec860 100644
--- a/roofit/roofit/src/RooNovosibirsk.cxx
+++ b/roofit/roofit/src/RooNovosibirsk.cxx
@@ -25,15 +25,15 @@ RooNovosibirsk implements the Novosibirsk function
 Function taken from H. Ikeda et al. NIM A441 (2000), p. 401 (Belle Collaboration)
 
 **/
-
+#include "RooNovosibirsk.h"
 #include "RooFit.h"
+#include "RooRealVar.h"
+#include "BatchHelpers.h"
+#include "RooVDTHeaders.h"
 
-#include <math.h>
 #include "TMath.h"
 
-#include "RooNovosibirsk.h"
-#include "RooRealVar.h"
-
+#include <cmath>
 using namespace std;
 
 ClassImp(RooNovosibirsk);
@@ -47,8 +47,8 @@ RooNovosibirsk::RooNovosibirsk(const char *name, const char *title,
   // parameter, respectively, as declared in the rdl file
   RooAbsPdf(name, title),
   x("x","x",this,_x),
-  width("width","width",this,_width),
   peak("peak","peak",this,_peak),
+  width("width","width",this,_width),
   tail("tail","tail",this,_tail)
 {
 }
@@ -58,8 +58,8 @@ RooNovosibirsk::RooNovosibirsk(const char *name, const char *title,
 RooNovosibirsk::RooNovosibirsk(const RooNovosibirsk& other, const char *name):
   RooAbsPdf(other,name),
   x("x",this,other.x),
-  width("width",this,other.width),
   peak("peak",this,other.peak),
+  width("width",this,other.width),
   tail("tail",this,other.tail)
 {
 }
@@ -92,6 +92,69 @@ Double_t RooNovosibirsk::evaluate() const
 
 ////////////////////////////////////////////////////////////////////////////////
 
+namespace {
+//Author: Emmanouil Michalainas, CERN 10 September 2019
+
+/* TMath::ASinH(x) needs to be replaced with ln( x + sqrt(x^2+1))
+ * argasinh -> the argument of TMath::ASinH()
+ * argln -> the argument of the logarithm that replaces AsinH
+ * asinh -> the value that the function evaluates to
+ * 
+ * ln is the logarithm that was solely present in the initial
+ * formula, that is before the asinh replacement
+ */
+template<class Tx, class Twidth, class Tpeak, class Ttail>
+void compute(	size_t batchSize, double * __restrict output,
+              Tx X, Tpeak P, Twidth W, Ttail T)
+{
+  constexpr double xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )
+  for (size_t i=0; i<batchSize; i++) {
+    double argasinh = 0.5*xi*T[i];
+    double argln = argasinh + 1/_rf_fast_isqrt(argasinh*argasinh +1);
+    double asinh = _rf_fast_log(argln);
+    
+    double argln2 = 1 -(X[i]-P[i])*T[i]/W[i];
+    double ln    = _rf_fast_log(argln2);
+    output[i] = ln/asinh;
+    output[i] *= -0.125*xi*xi*output[i];
+    output[i] -= 2.0/xi/xi*asinh*asinh;
+  }
+  
+  //faster if you exponentiate in a seperate loop (dark magic!)
+  for (size_t i=0; i<batchSize; i++) {
+    output[i] = _rf_fast_exp(output[i]);
+  }
+}
+};
+
+RooSpan<double> RooNovosibirsk::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
+  using namespace BatchHelpers;
+
+  EvaluateInfo info = getInfo( {&x, &peak, &width, &tail}, begin, batchSize );
+  if (info.nBatches == 0) {
+    return {};
+  }
+  auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
+  auto xData = x.getValBatch(begin, info.size);
+
+  if (info.nBatches==1 && !xData.empty()) {
+    compute(info.size, output.data(), xData.data(),
+    BracketAdapter<double> (peak),
+    BracketAdapter<double> (width),
+    BracketAdapter<double> (tail));
+  }
+  else {
+    compute(info.size, output.data(),
+    BracketAdapterWithMask (x,x.getValBatch(begin,info.size)),
+    BracketAdapterWithMask (peak,peak.getValBatch(begin,info.size)),
+    BracketAdapterWithMask (width,width.getValBatch(begin,info.size)),
+    BracketAdapterWithMask (tail,tail.getValBatch(begin,info.size)));
+  }
+  return output;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
 Int_t RooNovosibirsk::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* ) const
 {
   if (matchArgs(allVars,analVars,x)) return 1 ;
diff --git a/roofit/roofit/src/RooVoigtian.cxx b/roofit/roofit/src/RooVoigtian.cxx
index 90739a12a25a14..4c8599fa2ce997 100644
--- a/roofit/roofit/src/RooVoigtian.cxx
+++ b/roofit/roofit/src/RooVoigtian.cxx
@@ -25,18 +25,16 @@ algorithm. Select the faster algorithm either in the constructor, or with
 the selectFastAlgorithm() method.
 **/
 
-#include <cmath>
-#include <complex>
-
-#include "RooFit.h"
-
-#include "Riostream.h"
-
 #include "RooVoigtian.h"
+#include "RooFit.h"
 #include "RooAbsReal.h"
 #include "RooRealVar.h"
 #include "RooMath.h"
+#include "BatchHelpers.h"
+#include "RooVDTHeaders.h"
 
+#include <cmath>
+#include <complex>
 using namespace std;
 
 ClassImp(RooVoigtian);
@@ -54,7 +52,7 @@ RooVoigtian::RooVoigtian(const char *name, const char *title,
   sigma("sigma","Gauss Width",this,_sigma),
   _doFast(doFast)
 {
-  _invRootPi= 1./sqrt(atan2(0.,-1.));
+
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -64,7 +62,7 @@ RooVoigtian::RooVoigtian(const RooVoigtian& other, const char* name) :
   width("width",this,other.width),sigma("sigma",this,other.sigma),
   _doFast(other._doFast)
 {
-  _invRootPi= 1./sqrt(atan2(0.,-1.));
+
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -98,6 +96,66 @@ Double_t RooVoigtian::evaluate() const
   } else {
     v = RooMath::faddeeva(z);
   }
-  return c*_invRootPi*v.real();
+  return c * v.real();
+}
 
+////////////////////////////////////////////////////////////////////////////////
+
+namespace {
+//Author: Emmanouil Michalainas, CERN 11 September 2019
+
+template<class Tx, class Tmean, class Twidth, class Tsigma>
+void compute(	size_t batchSize, double * __restrict output,
+              Tx X, Tmean M, Twidth W, Tsigma S)
+{
+  constexpr double invSqrt2 = 0.707106781186547524400844362105;
+  for (size_t i=0; i<batchSize; i++) {
+    const double arg = (X[i]-M[i])*(X[i]-M[i]);
+    if (S[i]==0.0 && W[i]==0.0) {
+      output[i] = 1.0;
+    } else if (S[i]==0.0) {
+      output[i] = 1/(arg+0.25*W[i]*W[i]);
+    } else if (W[i]==0.0) {
+      output[i] = _rf_fast_exp(-0.5*arg/(S[i]*S[i]));
+    } else {
+      output[i] = invSqrt2/S[i];
+    }
+  }
+  
+  for (size_t i=0; i<batchSize; i++) {
+    if (S[i]!=0.0 && W[i]!=0.0) {
+      if (output[i] < 0) output[i] = -output[i];
+      const double factor = W[i]>0.0 ? 0.5 : -0.5;
+      std::complex<Double_t> z( output[i]*(X[i]-M[i]) , factor*output[i]*W[i] );
+      output[i] *= RooMath::faddeeva(z).real();
+    }
+  }
+}
+};
+
+RooSpan<double> RooVoigtian::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
+  using namespace BatchHelpers;
+
+  EvaluateInfo info = getInfo( {&x, &mean, &width, &sigma}, begin, batchSize );
+  if (info.nBatches == 0) {
+    return {};
+  }
+  auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
+  auto xData = x.getValBatch(begin, info.size);
+
+  if (info.nBatches==1 && !xData.empty()) {
+    compute(info.size, output.data(), xData.data(),
+    BracketAdapter<double> (mean),
+    BracketAdapter<double> (width),
+    BracketAdapter<double> (sigma));
+  }
+  else {
+    compute(info.size, output.data(),
+    BracketAdapterWithMask (x,x.getValBatch(begin,info.size)),
+    BracketAdapterWithMask (mean,mean.getValBatch(begin,info.size)),
+    BracketAdapterWithMask (width,width.getValBatch(begin,info.size)),
+    BracketAdapterWithMask (sigma,sigma.getValBatch(begin,info.size)));
+  }
+  return output;
 }
+
diff --git a/roofit/roofitcore/inc/RooVDTHeaders.h b/roofit/roofitcore/inc/RooVDTHeaders.h
index 39ec107abb8866..b0456469b6ed41 100644
--- a/roofit/roofitcore/inc/RooVDTHeaders.h
+++ b/roofit/roofitcore/inc/RooVDTHeaders.h
@@ -25,6 +25,7 @@
 #if defined(R__HAS_VDT)
 #include "vdt/exp.h"
 #include "vdt/log.h"
+#include "vdt/sqrt.h"
 
 inline double _rf_fast_exp(double x) {
   return vdt::fast_exp(x);
@@ -34,6 +35,10 @@ inline double _rf_fast_log(double x) {
   return vdt::fast_log(x);
 }
 
+inline double _rf_fast_isqrt(double x) {
+  return vdt::fast_isqrt(x);
+}
+
 #else
 #include <cmath>
 
@@ -45,6 +50,10 @@ inline double _rf_fast_log(double x) {
   return std::exp(x);
 }
 
+inline double _rf_fast_isqrt(double x) {
+  return 1/std::sqrt(x);
+}
+
 #endif