[RF] Vectorized more PDFs.

- DstD0BG
- Lognormal
- Novosibirsk
- Uniform
- batch evaluation of Voigtian
- Added isqrt in VDTHeaders

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:
 

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

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

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:
 

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

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 ;

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

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 ;