Logo ROOT   6.14/05
Reference Guide
TF1Convolution.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Authors: Lorenzo Moneta, AurĂ©lie Flandi 27/08/14
3 //
4 /**********************************************************************
5  * *
6  * Copyright (c) 2015 ROOT Team, CERN/PH-SFT *
7  * *
8  * *
9  **********************************************************************/
10 
11 #include "TF1Convolution.h"
12 #include "Riostream.h"
13 #include "TROOT.h"
14 #include "TObject.h"
15 #include "TObjString.h"
16 #include "TMath.h"
17 #include "Math/Integrator.h"
19 #include "Math/IntegratorOptions.h"
20 #include "Math/GaussIntegrator.h"
23 #include "Math/Functor.h"
24 #include "TVirtualFFT.h"
25 #include "TClass.h"
26 
27 /** \class TF1Convolution
28  \ingroup Hist
29  \brief Class wrapping convolution of two functions
30 
31 Class wrapping convolution of two functions: evaluation of \f$\int f(x)g(x-t)dx\f$
32 
33 The convolution is performed by default using FFTW if it is available .
34 One can pass optionally the range of the convolution (by default the first function range is used).
35 Note that when using Discrete Fourier Transform (as FFTW), it is a circular transform, so the functions should be
36 approximately zero at the end of the range. If they are significantly different than zero on one side (e.g. the left side)
37 a spill over will occur on the other side (e.g right side).
38 If no function range is given by default the function1 range + 10% is used
39 One should use also a not too small number of points for the DFT (a minimum of 1000). By default 10000 points are used.
40 */
41 
43 
44 class TF1Convolution_EvalWrapper
45 {
46  std::unique_ptr<TF1> fFunction1;
47  std::unique_ptr<TF1> fFunction2;
48  Double_t fT0;
49 
50 public:
51  TF1Convolution_EvalWrapper(std::unique_ptr<TF1> &f1, std::unique_ptr<TF1> &f2, Double_t t) : fT0(t)
52  {
53  fFunction1 = std::unique_ptr<TF1>((TF1 *)f1->Clone());
54  fFunction2 = std::unique_ptr<TF1>((TF1 *)f2->Clone());
55  }
57  {
58  // use EvalPar that is faster
59  Double_t xx[2];
60  xx[0] = x;
61  xx[1] = fT0-x;
62  return fFunction1->EvalPar(xx,nullptr) * fFunction2->EvalPar(xx+1,nullptr);
63  }
64 };
65 
66 ////////////////////////////////////////////////////////////////////////////////
67 /// Use copy instead of Clone
68 
69 void TF1Convolution::InitializeDataMembers(TF1* function1, TF1* function2, Bool_t useFFT)
70 {
71  if (function1) {
72  TF1 * fnew1 = (TF1*) function1->IsA()->New();
73  function1->Copy(*fnew1);
74  fFunction1 = std::unique_ptr<TF1>(fnew1);
75  }
76  if (function2) {
77  TF1 * fnew2 = (TF1*) function2->IsA()->New();
78  function2->Copy(*fnew2);
79  fFunction2 = std::unique_ptr<TF1>(fnew2);
80  }
81  if (fFunction1.get() == nullptr|| fFunction2.get() == nullptr)
82  Fatal("InitializeDataMembers","Invalid functions - Abort");
83 
84  // Set kNotGlobal bit
85  fFunction1->SetBit(TF1::kNotGlobal, kTRUE);
86  fFunction2->SetBit(TF1::kNotGlobal, kTRUE);
87 
88  // add by default an extra 10% on each side
89  fFunction1->GetRange(fXmin, fXmax);
90  Double_t range = fXmax - fXmin;
91  fXmin -= 0.1*range;
92  fXmax += 0.1*range;
93  fNofParams1 = fFunction1->GetNpar();
94  fNofParams2 = fFunction2->GetNpar();
95  fParams1 = std::vector<Double_t>(fNofParams1);
96  fParams2 = std::vector<Double_t>(fNofParams2);
97  fCstIndex = (fFunction1->GetParNumber("Constant") == -1)
98  ? -1
99  : fFunction2->GetParNumber("Constant"); // TODO: add dropConstantParam flag?
100  fFlagFFT = useFFT;
101  fFlagGraph = false;
102  fNofPoints = 10000;
103 
104  fParNames.reserve( fNofParams1 + fNofParams2);
105  for (int i=0; i<fNofParams1; i++)
106  {
107  fParams1[i] = fFunction1 -> GetParameter(i);
108  fParNames.push_back(fFunction1 -> GetParName(i) );
109  }
110  for (int i=0; i<fNofParams2; i++)
111  {
112  fParams2[i] = fFunction2 -> GetParameter(i);
113  if (i != fCstIndex) fParNames.push_back(fFunction2 -> GetParName(i) );
114  }
115  if (fCstIndex!=-1)
116  {
117  fFunction2 -> FixParameter(fCstIndex,1.);
118  fNofParams2 = fNofParams2-1;
119  fParams2.erase(fParams2.begin()+fCstIndex);
120  }
121 }
122 
123 ////////////////////////////////////////////////////////////////////////////////
124 /// constructor without arguments
125 
127 {
128  // Nothing to do
129 }
130 
131 ////////////////////////////////////////////////////////////////////////////////
132 /// constructor from the two function pointer and a flag is using FFT
133 
134 TF1Convolution::TF1Convolution(TF1* function1, TF1* function2, Bool_t useFFT)
135 {
136  InitializeDataMembers(function1,function2, useFFT);
137 }
138 
139 ////////////////////////////////////////////////////////////////////////////////
140 /// Constructor from the two function pointer and the convolution range
141 
143 {
144  InitializeDataMembers(function1, function2,useFFT);
145  if (xmin < xmax) {
146  fXmin = xmin;
147  fXmax = xmax;
148  } else {
149  Info("TF1Convolution", "Using default range [-inf, inf] for TF1Convolution");
150  SetRange(-TMath::Infinity(), TMath::Infinity());
151  }
152 }
153 
154 ////////////////////////////////////////////////////////////////////////////////
155 /// Constructor from a formula expression as f1 * f2 where f1 and f2 are two functions known to ROOT
156 
158 {
160 
161  TObjArray *objarray = formula.Tokenize("*");
162  std::vector < TString > stringarray(2);
163  std::vector < TF1* > funcarray(2);
164  for (int i=0; i<2; i++)
165  {
166  stringarray[i] = ((TObjString*)((*objarray)[i])) -> GetString();
167  stringarray[i].ReplaceAll(" ","");
168  funcarray[i] = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(stringarray[i]));
169  // case function is not found try to use as a TFormula
170  if (funcarray[i] == nullptr) {
171  TF1 * f = new TF1(TString::Format("f_conv_%d",i+1),stringarray[i]);
172  if (!f->GetFormula()->IsValid() )
173  Error("TF1Convolution","Invalid formula : %s",stringarray[i].Data() );
174  if (i == 0)
175  fFunction1 = std::unique_ptr<TF1>(f);
176  else
177  fFunction2 = std::unique_ptr<TF1>(f);
178  }
179  }
180  InitializeDataMembers(funcarray[0], funcarray[1],useFFT);
181  if (xmin < xmax) {
182  fXmin = xmin;
183  fXmax = xmax;
184  } else {
185  Info("TF1Convolution", "Using default range [-inf, inf] for TF1Convolution");
186  SetRange(-TMath::Infinity(), TMath::Infinity());
187  }
188 }
189 
190 ////////////////////////////////////////////////////////////////////////////////
191 /// constructor from 2 function names where f1 and f2 are two functions known to
192 /// ROOT
193 ///
194 /// if the function names are not known to ROOT, tries to interpret them as
195 /// TFormula
197 {
199  (TString)formula1.ReplaceAll(" ","");
200  (TString)formula2.ReplaceAll(" ","");
201  TF1* f1 = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(formula1));
202  TF1* f2 = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(formula2));
203  // if function do not exists try using TFormula
204  if (!f1) {
205  fFunction1 = std::unique_ptr<TF1>(new TF1("f_conv_1", formula1));
206  if (!fFunction1->GetFormula()->IsValid() )
207  Error("TF1Convolution","Invalid formula for : %s",formula1.Data() );
208  }
209  if (!f2) {
210  fFunction2 = std::unique_ptr<TF1>(new TF1("f_conv_1", formula2));
211  if (!fFunction2->GetFormula()->IsValid() )
212  Error("TF1Convolution","Invalid formula for : %s",formula2.Data() );
213  }
214  // if f1 or f2 are null ptr are not used in InitializeDataMembers
215  InitializeDataMembers(f1, f2,useFFT);
216  if (xmin < xmax) {
217  fXmin = xmin;
218  fXmax = xmax;
219  } else {
220  Info("TF1Convolution", "Using default range [-inf, inf] for TF1Convolution");
221  SetRange(-TMath::Infinity(), TMath::Infinity());
222  }
223 }
224 
225 ////////////////////////////////////////////////////////////////////////////////
226 /// Copy constructor (necessary to hold unique_ptr as member variable)
227 
229 {
230  conv.Copy((TObject &)*this);
231 }
232 
233 ////////////////////////////////////////////////////////////////////////////////
234 /// Operator =
235 
237 {
238  if (this != &rhs)
239  rhs.Copy(*this);
240  return *this;
241 }
242 
243 ////////////////////////////////////////////////////////////////////////////////
244 /// Perform the FFT of the two functions
245 
247 {
248  if (gDebug)
249  Info("MakeFFTConv","Making FFT convolution using %d points in range [%g,%g]",fNofPoints,fXmin,fXmax);
250 
251  std::vector < Double_t > x (fNofPoints);
252  std::vector < Double_t > in1(fNofPoints);
253  std::vector < Double_t > in2(fNofPoints);
254 
255  TVirtualFFT *fft1 = TVirtualFFT::FFT(1, &fNofPoints, "R2C K");
256  TVirtualFFT *fft2 = TVirtualFFT::FFT(1, &fNofPoints, "R2C K");
257  if (fft1 == nullptr || fft2 == nullptr) {
258  Warning("MakeFFTConv","Cannot use FFT, probably FFTW package is not available. Switch to numerical convolution");
259  fFlagFFT = false;
260  return;
261  }
262 
263  // apply a shift in order to have the second function centered around middle of the range of the convolution
264  Double_t shift2 = 0.5*(fXmin+fXmax);
265  Double_t x2;
266  for (int i=0; i<fNofPoints; i++)
267  {
268  x[i] = fXmin + (fXmax-fXmin)/(fNofPoints-1)*i;
269  x2 = x[i] - shift2;
270  in1[i] = fFunction1 -> EvalPar( &x[i], nullptr);
271  in2[i] = fFunction2 -> EvalPar( &x2, nullptr);
272  fft1 -> SetPoint(i, in1[i]);
273  fft2 -> SetPoint(i, in2[i]);
274  }
275  fft1 -> Transform();
276  fft2 -> Transform();
277 
278  //inverse transformation of the product
279 
280  TVirtualFFT *fftinverse = TVirtualFFT::FFT(1, &fNofPoints, "C2R K");
281  Double_t re1, re2, im1, im2, out_re, out_im;
282 
283  for (int i=0;i<=fNofPoints/2.;i++)
284  {
285  fft1 -> GetPointComplex(i,re1,im1);
286  fft2 -> GetPointComplex(i,re2,im2);
287  out_re = re1*re2 - im1*im2;
288  out_im = re1*im2 + re2*im1;
289  fftinverse -> SetPoint(i, out_re, out_im);
290  }
291  fftinverse -> Transform();
292 
293  // fill a graph with the result of the convolution
294  if (!fGraphConv)
295  fGraphConv = std::unique_ptr<TGraph>(new TGraph(fNofPoints));
296 
297  for (int i=0;i<fNofPoints;i++)
298  {
299  // we need this since we have applied a shift in the middle of f2
300  int j = i + fNofPoints/2;
301  if (j >= fNofPoints) j -= fNofPoints;
302  // need to normalize by dividing by the number of points and multiply by the bin width = Range/Number of points
303  fGraphConv->SetPoint(i, x[i], fftinverse->GetPointReal(j)*(fXmax-fXmin)/(fNofPoints*fNofPoints) );
304  }
305  fGraphConv->SetBit(TGraph::kIsSortedX); // indicate that points are sorted in X to speed up TGraph::Eval
306  fFlagGraph = true; // we can use the graph
307 
308  // delete the fft objects
309  delete fft1;
310  delete fft2;
311  delete fftinverse;
312 }
313 
314 ////////////////////////////////////////////////////////////////////////////////
315 
317 {
318  if (!fFlagGraph) MakeFFTConv();
319  // if cannot make FFT use numconv
320  if (fGraphConv)
321  return fGraphConv -> Eval(t);
322  else
323 
324  return EvalNumConv(t);
325 }
326 
327 ////////////////////////////////////////////////////////////////////////////////
328 /// Perform numerical convolution.
329 /// Could in principle cache the integral in a Graph as it is done for the FFTW
330 
332 {
333  TF1Convolution_EvalWrapper fconv( fFunction1, fFunction2, t);
334  Double_t result = 0;
335 
337  if (fXmin != - TMath::Infinity() && fXmax != TMath::Infinity() )
338  result = integrator.Integral(fXmin, fXmax);
339  else if (fXmin == - TMath::Infinity() && fXmax != TMath::Infinity() )
340  result = integrator.IntegralLow(fXmax);
341  else if (fXmin != - TMath::Infinity() && fXmax == TMath::Infinity() )
342  result = integrator.IntegralUp(fXmin);
343  else if (fXmin == - TMath::Infinity() && fXmax == TMath::Infinity() )
344  result = integrator.Integral();
345 
346  return result;
347 }
348 
349 ////////////////////////////////////////////////////////////////////////////////
350 /// Used in TF1 when doing the fit, will be evaluated at each point.
351 
353 {
354  if (p!=0) TF1Convolution::SetParameters(p); // first refresh the parameters
355 
356  Double_t result = 0.;
357  if (fFlagFFT)
358  result = EvalFFTConv(x[0]);
359  else
360  result = EvalNumConv(x[0]);
361  return result;
362 }
363 
364 ////////////////////////////////////////////////////////////////////////////////
365 
367 {
368  if (n<0) return;
369  fNofPoints = n;
370  if (fGraphConv) fGraphConv -> Set(fNofPoints);
371  fFlagGraph = false; // to indicate we need to re-do the graph
372 }
373 
374 ////////////////////////////////////////////////////////////////////////////////
375 
377 {
378  bool equalParams = true;
379  for (int i=0; i<fNofParams1; i++) {
380  fFunction1->SetParameter(i, params[i]);
381  equalParams &= (fParams1[i] == params[i]);
382  fParams1[i] = params[i];
383  }
384  Int_t k = 0;
385  Int_t offset = 0;
386  Int_t offset2 = 0;
387  if (fCstIndex!=-1) offset = 1;
388  Int_t totalnofparams = fNofParams1+fNofParams2+offset;
389  for (int i=fNofParams1; i<totalnofparams; i++) {
390  if (k==fCstIndex)
391  {
392  k++;
393  offset2=1;
394  continue;
395  }
396  fFunction2->SetParameter(k, params[i - offset2]);
397  equalParams &= (fParams2[k - offset2] == params[i - offset2]);
398  fParams2[k - offset2] = params[i - offset2];
399  k++;
400  }
401 
402  if (!equalParams) fFlagGraph = false; // to indicate we need to re-do the convolution
403 }
404 
405 ////////////////////////////////////////////////////////////////////////////////
406 
408  Double_t p4, Double_t p5, Double_t p6, Double_t p7)
409 {
410  Double_t params[]={p0,p1,p2,p3,p4,p5,p6,p7};
412 }
413 
414 ////////////////////////////////////////////////////////////////////////////////
415 
417 {
418  if (percentage<0) return;
419  double range = fXmax = fXmin;
420  fXmin -= percentage * range;
421  fXmax += percentage * range;
422  fFlagGraph = false; // to indicate we need to re-do the convolution
423 }
424 
425 ////////////////////////////////////////////////////////////////////////////////
426 
428 {
429  if (a >= b) {
430  Warning("SetRange", "Invalid range: %f >= %f", a, b);
431  return;
432  }
433 
434  fXmin = a;
435  fXmax = b;
436  if (fFlagFFT && ( a==-TMath::Infinity() || b==TMath::Infinity() ) )
437  {
438  Warning("TF1Convolution::SetRange()","In FFT mode, range can not be infinite. Infinity has been replaced by range of first function plus a bufferzone to avoid spillover.");
439  if (a ==-TMath::Infinity()) fXmin = fFunction1 -> GetXmin();
440  if ( b== TMath::Infinity()) fXmax = fFunction1 -> GetXmax();
441  // add a spill over of 10% in this case
442  SetExtraRange(0.1);
443  }
444  fFlagGraph = false; // to indicate we need to re-do the convolution
445 }
446 
447 ////////////////////////////////////////////////////////////////////////////////
448 
450 {
451  a = fXmin;
452  b = fXmax;
453 }
454 
455 ////////////////////////////////////////////////////////////////////////////////
456 /// Update the two component functions of the convolution
457 
459 {
460  fFunction1->Update();
461  fFunction2->Update();
462 }
463 
464 ////////////////////////////////////////////////////////////////////////////////
465 
467 {
468  // copy numbers
469  ((TF1Convolution &)obj).fXmin = fXmin;
470  ((TF1Convolution &)obj).fXmax = fXmax;
471  ((TF1Convolution &)obj).fNofParams1 = fNofParams1;
472  ((TF1Convolution &)obj).fNofParams2 = fNofParams2;
473  ((TF1Convolution &)obj).fCstIndex = fCstIndex;
474  ((TF1Convolution &)obj).fNofPoints = fNofPoints;
475  ((TF1Convolution &)obj).fFlagFFT = fFlagFFT;
476  ((TF1Convolution &)obj).fFlagGraph = false; // since we're not copying the graph
477 
478  // copy vectors
479  ((TF1Convolution &)obj).fParams1 = fParams1;
480  ((TF1Convolution &)obj).fParams2 = fParams2;
481  ((TF1Convolution &)obj).fParNames = fParNames;
482 
483  // Clone unique_ptr's
484  ((TF1Convolution &)obj).fFunction1 = std::unique_ptr<TF1>((TF1 *)fFunction1->Clone());
485  ((TF1Convolution &)obj).fFunction2 = std::unique_ptr<TF1>((TF1 *)fFunction2->Clone());
486  // fGraphConv is transient anyway, so we don't bother to copy it
487 }
An array of TObjects.
Definition: TObjArray.h:37
float xmin
Definition: THbookFile.cxx:93
virtual TFormula * GetFormula()
Definition: TF1.h:437
static double p3(double t, double a, double b, double c, double d)
virtual void Copy(TObject &f1) const
Copy this F1 to a new F1.
Definition: TF1.cxx:906
void Update()
Update the two component functions of the convolution.
void Fatal(const char *location, const char *msgfmt,...)
Collectable string class.
Definition: TObjString.h:28
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
TF1Convolution()
constructor without arguments
#define gROOT
Definition: TROOT.h:410
Basic string class.
Definition: TString.h:131
#define f(i)
Definition: RSha256.hxx:104
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Double_t operator()(const Double_t *x, const Double_t *p)
Used in TF1 when doing the fit, will be evaluated at each point.
TF1Convolution & operator=(const TF1Convolution &rhs)
Operator =.
TRObject operator()(const T1 &t1) const
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
Class wrapping convolution of two functions.
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
Definition: TString.cxx:2286
static double p2(double t, double a, double b, double c)
void MakeFFTConv()
Perform the FFT of the two functions.
void Info(const char *location, const char *msgfmt,...)
void InitializeDataMembers(TF1 *function1, TF1 *function2, Bool_t useFFT)
Use copy instead of Clone.
static IntegrationOneDim::Type DefaultIntegratorType()
Double_t Infinity()
Returns an infinity as defined by the IEEE standard.
Definition: TMath.h:913
void Error(const char *location, const char *msgfmt,...)
double IntegralLow(const IGenFunction &f, double b)
evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b)
Definition: Integrator.h:292
double Integral(Function &f, double a, double b)
evaluate the Integral of a function f over the defined interval (a,b)
Definition: Integrator.h:496
void SetRange(Double_t a, Double_t b)
double IntegralUp(const IGenFunction &f, double a)
evaluate the Integral of a function f over the semi-infinite interval (a,+inf)
Definition: Integrator.h:274
virtual Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const =0
auto * a
Definition: textangle.C:12
static TVirtualFFT * FFT(Int_t ndim, Int_t *n, Option_t *option)
Returns a pointer to the FFT of requested size and type.
Double_t EvalNumConv(Double_t t)
Perform numerical convolution.
void SetNofPointsFFT(Int_t n)
User Class for performing numerical integration of a function in one dimension.
Definition: Integrator.h:94
static void InitStandardFunctions()
Create the basic function objects.
Definition: TF1.cxx:2376
static double p1(double t, double a, double b)
float xmax
Definition: THbookFile.cxx:93
virtual void Update()
Called by functions such as SetRange, SetNpx, SetParameters to force the deletion of the associated h...
Definition: TF1.cxx:3575
void Warning(const char *location, const char *msgfmt,...)
TVirtualFFT is an interface class for Fast Fourier Transforms.
Definition: TVirtualFFT.h:88
TObjArray * Tokenize(const TString &delim) const
This function is used to isolate sequential tokens in a TString.
Definition: TString.cxx:2172
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
graph is sorted in X points
Definition: TGraph.h:73
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t EvalFFTConv(Double_t t)
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TNamed.cxx:74
Mother of all ROOT objects.
Definition: TObject.h:37
virtual Int_t GetNpar() const
Definition: TF1.h:465
void Copy(TObject &obj) const
Copy this to obj.
virtual void GetRange(Double_t *xmin, Double_t *xmax) const
Return range of a generic N-D function.
Definition: TF1.cxx:2170
1-Dim function class
Definition: TF1.h:211
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
TF1 * f1
Definition: legend1.C:11
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
R__EXTERN Int_t gDebug
Definition: Rtypes.h:86
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:618
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=0)
Evaluate function with given coordinates and parameters.
Definition: TF1.cxx:1365
Bool_t IsValid() const
Definition: TFormula.h:204
const Bool_t kTRUE
Definition: RtypesCore.h:87
const Int_t n
Definition: legend1.C:16
void GetRange(Double_t &a, Double_t &b) const
virtual Int_t GetParNumber(const char *name) const
Definition: TF1.h:517
void SetExtraRange(Double_t percentage)
const char * Data() const
Definition: TString.h:364
void SetParameters(const Double_t *params)