diff --git a/hist/hist/inc/TF1.h b/hist/hist/inc/TF1.h
index 7bfa97f27bc8aaae9d095fdcf45236f180f3f4bc..0d2e6abbac24041925a7cbb4ccb9ad6b7d937895 100644
--- a/hist/hist/inc/TF1.h
+++ b/hist/hist/inc/TF1.h
@@ -24,6 +24,7 @@
#include "RConfigure.h"
#include <functional>
#include <cassert>
+#include <memory> // for shared_ptr
#include "TFormula.h"
#include "TAttLine.h"
#include "TAttFill.h"
@@ -37,6 +38,7 @@ class TF1;
class TH1;
class TAxis;
class TMethodCall;
+class TF1NormSum;
namespace ROOT {
namespace Fit {
@@ -257,14 +259,22 @@ protected:
TF1Parameters *fParams = nullptr; //Pointer to Function parameters object (exists only for not-formula functions)
/// General constructor for TF1. Most of the other constructors delegate on it
- TF1(EFType functionType, const char *name, Double_t xmin, Double_t xmax, Int_t npar, Int_t ndim, EAddToList addToGlobList, TF1Parameters *params = nullptr, TF1FunctorPointer * functor = nullptr):
- TNamed(name, name), TAttLine(), TAttFill(), TAttMarker(), fXmin(xmin), fXmax(xmax), fNpar(npar), fNdim(ndim),
+ TF1(EFType functionType, const char *name, Double_t xmin, Double_t xmax, Int_t npar, Int_t ndim, EAddToList addToGlobList, TF1Parameters *params = nullptr, TF1FunctorPointer * functor = nullptr):
+ TNamed(name, name), TAttLine(), TAttFill(), TAttMarker(), fXmin(xmin), fXmax(xmax), fNpar(npar), fNdim(ndim),
fType(functionType), fParErrors(npar), fParMin(npar), fParMax(npar), fFunctor(functor), fParams(params)
{
DoInitialize(addToGlobList);
};
+private:
+ // NSUM parsing helper functions
+ void DefineNSUMTerm(TObjArray *newFuncs, TObjArray *coeffNames,
+ TString &fullFormula,
+ TString &formula, int termStart, int termEnd,
+ Double_t xmin, Double_t xmax);
+ int TermCoeffLength(TString &term);
+
public:
template <class T>
@@ -305,7 +315,6 @@ public:
TF1(const char *name, const char *formula, Double_t xmin = 0, Double_t xmax = 1, EAddToList addToGlobList = EAddToList::kDefault);
TF1(const char *name, Double_t xmin, Double_t xmax, Int_t npar, Int_t ndim = 1, EAddToList addToGlobList = EAddToList::kDefault);
TF1(const char *name, Double_t (*fcn)(Double_t *, Double_t *), Double_t xmin = 0, Double_t xmax = 1, Int_t npar = 0, Int_t ndim = 1, EAddToList addToGlobList = EAddToList::kDefault);
- TF1(const char *name, Double_t (*fcn)(const Double_t *, const Double_t *), Double_t xmin = 0, Double_t xmax = 1, Int_t npar = 0, Int_t ndim = 1, EAddToList addToGlobList = EAddToList::kDefault);
template <class T>
TF1(const char *name, std::function<T(const T *data, const Double_t *param)> &fcn, Double_t xmin = 0, Double_t xmax = 1, Int_t npar = 0, Int_t ndim = 1, EAddToList addToGlobList = EAddToList::kDefault):
@@ -316,17 +325,18 @@ public:
}
////////////////////////////////////////////////////////////////////////////////
- /// Constructor using a pointer to function.
+ /// Constructor using a pointer to real function.
///
/// \param npar is the number of free parameters used by the function
///
/// This constructor creates a function of type C when invoked
/// with the normal C++ compiler.
///
+ /// see test program test/stress.cxx (function stress1) for an example.
+ /// note the interface with an intermediate pointer.
///
/// WARNING! A function created with this constructor cannot be Cloned
-
template <class T>
TF1(const char *name, T(*fcn)(const T *, const Double_t *), Double_t xmin = 0, Double_t xmax = 1, Int_t npar = 0, Int_t ndim = 1, EAddToList addToGlobList = EAddToList::kDefault):
TF1(EFType::kTemplVec, name, xmin, xmax, npar, ndim, addToGlobList, new TF1Parameters(npar), new TF1FunctorPointerImpl<T>(fcn))
diff --git a/hist/hist/src/TF1.cxx b/hist/hist/src/TF1.cxx
index a1c00176df75e2d70f4591955662565a6f75b595..ffeb9e20a508f97631a7f4613b0c3b15faa199d0 100644
--- a/hist/hist/src/TF1.cxx
+++ b/hist/hist/src/TF1.cxx
@@ -25,6 +25,7 @@
#include "TClass.h"
#include "TMethodCall.h"
#include "TF1Helper.h"
+#include "TF1NormSum.h"
#include "TVirtualMutex.h"
#include "Math/WrappedFunction.h"
#include "Math/WrappedTF1.h"
@@ -431,9 +432,65 @@ TF1::TF1(const char *name, const char *formula, Double_t xmin, Double_t xmax, EA
fXmax = xmin;
}
// create rep formula (no need to add to gROOT list since we will add the TF1 object)
- fFormula = new TFormula(name, formula, false);
- fNpar = fFormula->GetNpar();
- fNdim = fFormula->GetNdim();
+ // First check if we need NSUM syntax:
+ if (TString(formula, 5) == "NSUM(" && formula[strlen(formula)-1] == ')') {
+ // using comma as delimiter
+ char delimiter = ',';
+ // first, remove "NSUM(" and ")" and spaces
+ TString formDense = TString(formula)(5,strlen(formula)-5-1);
+ formDense.ReplaceAll(' ', "");
+
+ // make sure standard functions are defined (e.g. gaus, expo)
+ InitStandardFunctions();
+
+ // Go char-by-char to split terms and define the relevant functions
+ int parenCount = 0;
+ int termStart = 0;
+ TObjArray *newFuncs = new TObjArray();
+ newFuncs->SetOwner(kTRUE);
+ TObjArray *coeffNames = new TObjArray();
+ coeffNames->SetOwner(kTRUE);
+ TString fullFormula("");
+ for (int i = 0; i < formDense.Length(); ++i) {
+ if (formDense[i] == '(')
+ parenCount++;
+ else if (formDense[i] == ')')
+ parenCount--;
+ else if (formDense[i] == delimiter && parenCount == 0) {
+ // term goes from termStart to i
+ DefineNSUMTerm(newFuncs, coeffNames, fullFormula, formDense, termStart, i, xmin, xmax);
+ termStart = i + 1;
+ }
+ }
+ DefineNSUMTerm(newFuncs, coeffNames, fullFormula, formDense, termStart, formDense.Length(), xmin, xmax);
+
+ TF1NormSum *normSum = new TF1NormSum(fullFormula, xmin, xmax);
+
+ fNpar = normSum->GetNpar();
+ fNdim = 1; // (note: may want to extend functionality in the future)
+
+ fType = EFType::kPtrScalarFreeFcn; // (note: may want to add new fType for this case)
+ // fFunctor = ROOT::Math::ParamFunctor(normSum);
+ using Fnc_t = typename ROOT::Internal::GetFunctorType<decltype(ROOT::Internal::GetTheRightOp(&TF1NormSum::operator()))>::type;
+ fFunctor = new TF1::TF1FunctorPointerImpl<Fnc_t>(ROOT::Math::ParamFunctorTempl<Fnc_t>(normSum));
+ fParams = new TF1Parameters(fNpar);
+ fParams->SetParameters(&(normSum->GetParameters())[0]); // inherit default parameters from normSum
+
+ // Parameter names
+ for (int i = 0; i < fNpar; i++) {
+ if (coeffNames->At(i) != nullptr) {
+ TString coeffName = ((TObjString *)coeffNames->At(i))->GetString();
+ this->SetParName(i, (const char *)coeffName);
+ } else {
+ this->SetParName(i, normSum->GetParName(i));
+ }
+ }
+
+ } else { // regular TFormula
+ fFormula = new TFormula(name, formula, false);
+ fNpar = fFormula->GetNpar();
+ fNdim = fFormula->GetNdim();
+ }
if (fNpar) {
fParErrors.resize(fNpar);
fParMin.resize(fNpar);
@@ -499,22 +556,6 @@ TF1::TF1(const char *name, Double_t (*fcn)(Double_t *, Double_t *), Double_t xmi
TF1(EFType::kPtrScalarFreeFcn, name, xmin, xmax, npar, ndim, addToGlobList, new TF1Parameters(npar), new TF1FunctorPointerImpl<double>(ROOT::Math::ParamFunctor(fcn)))
{}
-////////////////////////////////////////////////////////////////////////////////
-/// Constructor using a pointer to real function.
-///
-/// \param npar is the number of free parameters used by the function
-///
-/// This constructor creates a function of type C when invoked
-/// with the normal C++ compiler.
-///
-/// see test program test/stress.cxx (function stress1) for an example.
-/// note the interface with an intermediate pointer.
-///
-/// WARNING! A function created with this constructor cannot be Cloned.
-
-TF1::TF1(const char *name, Double_t (*fcn)(const Double_t *, const Double_t *), Double_t xmin, Double_t xmax, Int_t npar, Int_t ndim, EAddToList addToGlobList) :
- TF1(EFType::kPtrScalarFreeFcn, name, xmin, xmax, npar, ndim, addToGlobList, new TF1Parameters(npar), new TF1FunctorPointerImpl<double>(ROOT::Math::ParamFunctor(fcn)))
-{}
////////////////////////////////////////////////////////////////////////////////
/// Constructor using the Functor class.
@@ -615,6 +656,67 @@ Bool_t TF1::AddToGlobalList(Bool_t on)
return prevStatus;
}
+////////////////////////////////////////////////////////////////////////////////
+/// Helper functions for NSUM parsing
+
+// Defines the formula that a given term uses, if not already defined,
+// and appends "sanitized" formula to `fullFormula` string
+void TF1::DefineNSUMTerm(TObjArray *newFuncs, TObjArray *coeffNames,
+ TString &fullFormula,
+ TString &formula, int termStart, int termEnd,
+ Double_t xmin, Double_t xmax) {
+ TString originalTerm = formula(termStart, termEnd-termStart);
+ int coeffLength = TermCoeffLength(originalTerm);
+ if (coeffLength != -1)
+ termStart += coeffLength + 1;
+
+ // `originalFunc` is the real formula and `cleanedFunc` is the
+ // sanitized version that will not confuse the TF1NormSum
+ // constructor
+ TString originalFunc = formula(termStart, termEnd-termStart);
+ TString cleanedFunc = TString(formula(termStart, termEnd-termStart))
+ .ReplaceAll('+', "<plus>")
+ .ReplaceAll('*',"<times>");
+
+ // define function (if necessary)
+ if (!gROOT->GetListOfFunctions()->FindObject(cleanedFunc))
+ newFuncs->Add(new TF1(cleanedFunc, originalFunc, xmin, xmax));
+
+ // append sanitized term to `fullFormula`
+ if (fullFormula.Length() != 0)
+ fullFormula.Append('+');
+
+ // include numerical coefficient
+ if (coeffLength != -1 && originalTerm[0] != '[')
+ fullFormula.Append(originalTerm(0, coeffLength+1));
+
+ // add coefficient name
+ if (coeffLength != -1 && originalTerm[0] == '[')
+ coeffNames->Add(new TObjString(TString(originalTerm(1,coeffLength-2))));
+ else
+ coeffNames->Add(nullptr);
+
+ fullFormula.Append(cleanedFunc);
+}
+
+
+// Returns length of coeff at beginning of a given term, not counting the '*'
+// Returns -1 if no coeff found
+// Coeff can be either a number or parameter name
+int TF1::TermCoeffLength(TString &term) {
+ int firstAsterisk = term.First('*');
+ if (firstAsterisk == -1) // no asterisk found
+ return -1;
+
+ if (TString(term(0,firstAsterisk)).IsFloat())
+ return firstAsterisk;
+
+ if (term[0] == '[' && term[firstAsterisk-1] == ']'
+ && TString(term(1,firstAsterisk-2)).IsAlnum())
+ return firstAsterisk;
+
+ return -1;
+}
////////////////////////////////////////////////////////////////////////////////
/// Operator =
diff --git a/hist/hist/test/CMakeLists.txt b/hist/hist/test/CMakeLists.txt
index a851ac9ea005de314b6cef0715c44eb028ebd83c..dec3662c0bcebf3a6b6c454dbd8a9ba747192268 100644
--- a/hist/hist/test/CMakeLists.txt
+++ b/hist/hist/test/CMakeLists.txt
@@ -1,2 +1,5 @@
ROOT_ADD_GTEST(testTProfile2Poly test_tprofile2poly.cxx LIBRARIES Hist Matrix MathCore RIO)
ROOT_ADD_GTEST(testTHn THn.cxx LIBRARIES Hist Matrix MathCore RIO)
+
+# ROOT_ADD_GTEST(testTF1 test_tf1.cxx LIBRARIES Core Hist Matrix MathCore RIO)
+ROOT_ADD_GTEST(testTF1 test_tf1.cxx LIBRARIES Hist)
diff --git a/hist/hist/test/test_tf1.cxx b/hist/hist/test/test_tf1.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..312bd3ef5045db4885180c2f02d7c72ea289add9
--- /dev/null
+++ b/hist/hist/test/test_tf1.cxx
@@ -0,0 +1,71 @@
+#include "TF1.h"
+#include "TObjString.h"
+// #include "TObject.h"
+
+#include "gtest/gtest.h"
+
+#include <iostream>
+
+using namespace std;
+
+Float_t delta = 0.00000000001;
+
+void coeffNamesGeneric(TString *formula, TObjArray *coeffNames) {
+ TF1 *cn0 = new TF1("cn0", *formula, 0, 1);
+ ASSERT_EQ(cn0->GetNpar(), coeffNames->GetEntries());
+ for (int i = 0; i < coeffNames->GetEntries(); i++) {
+ TObjString *coeffObj = (TObjString *)coeffNames->At(i);
+ TString coeffName = coeffObj->GetString();
+ EXPECT_EQ(coeffName, TString(cn0->GetParName(i)));
+ }
+}
+
+void test_coeffNames() {
+ // cout << "About to start" << endl;
+
+ TObjArray *coeffNames = new TObjArray();
+ coeffNames->SetOwner(kTRUE);
+ coeffNames->Add(new TObjString("sg"));
+ coeffNames->Add(new TObjString("bg"));
+ coeffNames->Add(new TObjString("Mean"));
+ coeffNames->Add(new TObjString("Sigma"));
+ coeffNames->Add(new TObjString("Slope"));
+ TString *formula = new TString("NSUM([sg] * gaus, [bg] * expo)");
+
+ // cout << "Almost done" << endl;
+
+ coeffNamesGeneric(formula, coeffNames);
+}
+
+void test_normalization() {
+ double xmin = -5;
+ double xmax = 5;
+ TF1 *n0 = new TF1("n0", "NSUM(.5 * gaus, .5 * (x+[0])**2)", xmin, xmax);
+ EXPECT_NEAR(n0->Integral(xmin, xmax), 1, delta);
+ n0->SetParameter(4,1); // should not affect integral
+ EXPECT_NEAR(n0->Integral(xmin, xmax), 1, delta);
+ n0->SetParameter(0,0);
+ EXPECT_NEAR(n0->Integral(xmin, xmax), .5, delta);
+
+ TF1 *n1 = new TF1("n1", "NSUM([sg] * gaus, [bg] * (x+[0])**2)", xmin, xmax);
+ n1->SetParameter(0,.5);
+ n1->SetParameter(1,.5);
+ EXPECT_NEAR(n1->Integral(xmin, xmax), 1, delta);
+ n0->SetParameter(0,0);
+ EXPECT_NEAR(n0->Integral(xmin, xmax), .5, delta);
+
+ TF1 *n2 = new TF1("n2", "NSUM([sg] * gaus, -0.5 * (x+[0])**2)", xmin, xmax);
+ n2->SetParameter(0,.5);
+ EXPECT_NEAR(n2->GetParameter(1), -.5, delta);
+ EXPECT_NEAR(n2->Integral(xmin, xmax), 0, delta);
+ n2->SetParameter(0,0);
+ EXPECT_NEAR(n2->Integral(xmin, xmax), -.5, delta);
+}
+
+TEST(TF1, CoeffNames) {
+ test_coeffNames();
+}
+
+TEST(TF1, Normalization) {
+ test_normalization();
+}