ROOT   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"
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
31Class wrapping convolution of two functions: evaluation of \f$\int f(x)g(x-t)dx\f$
32
33The convolution is performed by default using FFTW if it is available .
34One can pass optionally the range of the convolution (by default the first function range is used).
35Note that when using Discrete Fourier Transform (as FFTW), it is a circular transform, so the functions should be
36approximately zero at the end of the range. If they are significantly different than zero on one side (e.g. the left side)
37a spill over will occur on the other side (e.g right side).
38If no function range is given by default the function1 range + 10% is used
39One 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
44class TF1Convolution_EvalWrapper
45{
46 TF1 * fFunc1;
47 TF1 * fFunc2;
48 Double_t fT0;
49
50public:
51 TF1Convolution_EvalWrapper(TF1 &f1, TF1 &f2, Double_t t) :
52 fFunc1(&f1),
53 fFunc2(&f2),
54 fT0(t)
55 {}
57 {
58 // use EvalPar that is faster
59 Double_t xx[2];
60 xx[0] = x;
61 xx[1] = fT0-x;
62 return fFunc1->EvalPar(xx,nullptr) * fFunc2->EvalPar(xx+1,nullptr);
63 }
64};
65
66////////////////////////////////////////////////////////////////////////////////
67/// Use copy instead of Clone
68
69void TF1Convolution::InitializeDataMembers(TF1* function1, TF1* function2, Bool_t useFFT)
70{
71 if (function1) {
72 // functions must be 1d- if not flag an error
73 if (function1->GetNdim() != 1)
74 Error("InitializeDataMembers","function1 %s is not of dimension 1 ",function1->GetName());
75 //TF1 * fnew1 = (TF1*) function1->IsA()->New();
76 // since function1 is a TF1 (cannot be a derived class) we can instantiate it directly
77 fFunction1 = std::unique_ptr<TF1> (new TF1());
78 function1->Copy(*fFunction1);
79 }
80 if (function2) {
81 if (function2->GetNdim() != 1)
82 Error("InitializeDataMembers","function2 %s is not of dimension 1 ",function2->GetName());
83 //TF1 * fnew2 = (TF1*) function2->IsA()->New();
84 fFunction2 = std::unique_ptr<TF1>(new TF1());
85 function2->Copy(*fFunction2);
86 }
87 if (fFunction1.get() == nullptr|| fFunction2.get() == nullptr)
88 Fatal("InitializeDataMembers","Invalid functions - Abort");
89
90 // Set kNotGlobal bit
93
94 // add by default an extra 10% on each side
95 fFunction1->GetRange(fXmin, fXmax);
96 Double_t range = fXmax - fXmin;
97 fXmin -= 0.1*range;
98 fXmax += 0.1*range;
99 fNofParams1 = fFunction1->GetNpar();
100 fNofParams2 = fFunction2->GetNpar();
101 fParams1 = std::vector<Double_t>(fNofParams1);
102 fParams2 = std::vector<Double_t>(fNofParams2);
103 fCstIndex = (fFunction1->GetParNumber("Constant") == -1)
104 ? -1
105 : fFunction2->GetParNumber("Constant"); // TODO: add dropConstantParam flag?
106 fFlagFFT = useFFT;
107 fFlagGraph = false;
108 fNofPoints = 10000;
109
111 for (int i=0; i<fNofParams1; i++)
112 {
113 fParams1[i] = fFunction1 -> GetParameter(i);
114 fParNames.push_back(fFunction1 -> GetParName(i) );
115 }
116 for (int i=0; i<fNofParams2; i++)
117 {
118 fParams2[i] = fFunction2 -> GetParameter(i);
119 if (i != fCstIndex) fParNames.push_back(fFunction2 -> GetParName(i) );
120 }
121 if (fCstIndex!=-1)
122 {
123 fFunction2 -> FixParameter(fCstIndex,1.);
125 fParams2.erase(fParams2.begin()+fCstIndex);
126 }
127}
128
129////////////////////////////////////////////////////////////////////////////////
130/// constructor without arguments
131
133{
134 // Nothing to do
135}
136
137////////////////////////////////////////////////////////////////////////////////
138/// constructor from the two function pointer and a flag is using FFT
139
140TF1Convolution::TF1Convolution(TF1* function1, TF1* function2, Bool_t useFFT)
141{
142 InitializeDataMembers(function1,function2, useFFT);
143}
144
145////////////////////////////////////////////////////////////////////////////////
146/// Constructor from the two function pointer and the convolution range
147
149{
150 InitializeDataMembers(function1, function2,useFFT);
151 if (xmin < xmax) {
152 fXmin = xmin;
153 fXmax = xmax;
154 } else {
155 Info("TF1Convolution", "Using default range [-inf, inf] for TF1Convolution");
157 }
158}
159
160////////////////////////////////////////////////////////////////////////////////
161/// Constructor from a formula expression as f1 * f2 where f1 and f2 are two functions known to ROOT
162
164{
166
167 TObjArray *objarray = formula.Tokenize("*");
168 std::vector < TString > stringarray(2);
169 std::vector < TF1* > funcarray(2);
170 for (int i=0; i<2; i++)
171 {
172 stringarray[i] = ((TObjString*)((*objarray)[i])) -> GetString();
173 stringarray[i].ReplaceAll(" ","");
174 funcarray[i] = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(stringarray[i]));
175 // case function is not found try to use as a TFormula
176 if (funcarray[i] == nullptr) {
177 TF1 * f = new TF1(TString::Format("f_conv_%d",i+1),stringarray[i]);
178 if (!f->GetFormula()->IsValid() )
179 Error("TF1Convolution","Invalid formula : %s",stringarray[i].Data() );
180 if (i == 0)
181 fFunction1 = std::unique_ptr<TF1>(f);
182 else
183 fFunction2 = std::unique_ptr<TF1>(f);
184 }
185 }
186 InitializeDataMembers(funcarray[0], funcarray[1],useFFT);
187 if (xmin < xmax) {
188 fXmin = xmin;
189 fXmax = xmax;
190 } else {
191 Info("TF1Convolution", "Using default range [-inf, inf] for TF1Convolution");
193 }
194}
195
196////////////////////////////////////////////////////////////////////////////////
197/// constructor from 2 function names where f1 and f2 are two functions known to
198/// ROOT
199///
200/// if the function names are not known to ROOT, tries to interpret them as
201/// TFormula
203{
205 (TString)formula1.ReplaceAll(" ","");
206 (TString)formula2.ReplaceAll(" ","");
207 TF1* f1 = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(formula1));
208 TF1* f2 = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(formula2));
209 // if function do not exists try using TFormula
210 if (!f1) {
211 fFunction1 = std::unique_ptr<TF1>(new TF1("f_conv_1", formula1));
212 if (!fFunction1->GetFormula()->IsValid() )
213 Error("TF1Convolution","Invalid formula for : %s",formula1.Data() );
214 }
215 if (!f2) {
216 fFunction2 = std::unique_ptr<TF1>(new TF1("f_conv_1", formula2));
217 if (!fFunction2->GetFormula()->IsValid() )
218 Error("TF1Convolution","Invalid formula for : %s",formula2.Data() );
219 }
220 // if f1 or f2 are null ptr are not used in InitializeDataMembers
221 InitializeDataMembers(f1, f2,useFFT);
222 if (xmin < xmax) {
223 fXmin = xmin;
224 fXmax = xmax;
225 } else {
226 Info("TF1Convolution", "Using default range [-inf, inf] for TF1Convolution");
228 }
229}
230
231////////////////////////////////////////////////////////////////////////////////
232/// Copy constructor (necessary to hold unique_ptr as member variable)
233
235{
236 conv.Copy((TObject &)*this);
237}
238
239////////////////////////////////////////////////////////////////////////////////
240/// Operator =
241
243{
244 if (this != &rhs)
245 rhs.Copy(*this);
246 return *this;
247}
248
249////////////////////////////////////////////////////////////////////////////////
250/// Perform the FFT of the two functions
251
253{
254 if (gDebug)
255 Info("MakeFFTConv","Making FFT convolution using %d points in range [%g,%g]",fNofPoints,fXmin,fXmax);
256
257 std::vector < Double_t > x (fNofPoints);
258 std::vector < Double_t > in1(fNofPoints);
259 std::vector < Double_t > in2(fNofPoints);
260
261 TVirtualFFT *fft1 = TVirtualFFT::FFT(1, &fNofPoints, "R2C K");
262 TVirtualFFT *fft2 = TVirtualFFT::FFT(1, &fNofPoints, "R2C K");
263 if (fft1 == nullptr || fft2 == nullptr) {
264 Warning("MakeFFTConv","Cannot use FFT, probably FFTW package is not available. Switch to numerical convolution");
265 fFlagFFT = false;
266 return;
267 }
268
269 // apply a shift in order to have the second function centered around middle of the range of the convolution
270 Double_t shift2 = 0.5*(fXmin+fXmax);
271 Double_t x2;
272 for (int i=0; i<fNofPoints; i++)
273 {
274 x[i] = fXmin + (fXmax-fXmin)/(fNofPoints-1)*i;
275 x2 = x[i] - shift2;
276 in1[i] = fFunction1 -> EvalPar( &x[i], nullptr);
277 in2[i] = fFunction2 -> EvalPar( &x2, nullptr);
278 fft1 -> SetPoint(i, in1[i]);
279 fft2 -> SetPoint(i, in2[i]);
280 }
281 fft1 -> Transform();
282 fft2 -> Transform();
283
284 //inverse transformation of the product
285
286 TVirtualFFT *fftinverse = TVirtualFFT::FFT(1, &fNofPoints, "C2R K");
287 Double_t re1, re2, im1, im2, out_re, out_im;
288
289 for (int i=0;i<=fNofPoints/2.;i++)
290 {
291 fft1 -> GetPointComplex(i,re1,im1);
292 fft2 -> GetPointComplex(i,re2,im2);
293 out_re = re1*re2 - im1*im2;
294 out_im = re1*im2 + re2*im1;
295 fftinverse -> SetPoint(i, out_re, out_im);
296 }
297 fftinverse -> Transform();
298
299 // fill a graph with the result of the convolution
300 if (!fGraphConv)
301 fGraphConv = std::unique_ptr<TGraph>(new TGraph(fNofPoints));
302
303 for (int i=0;i<fNofPoints;i++)
304 {
305 // we need this since we have applied a shift in the middle of f2
306 int j = i + fNofPoints/2;
307 if (j >= fNofPoints) j -= fNofPoints;
308 // need to normalize by dividing by the number of points and multiply by the bin width = Range/Number of points
309 fGraphConv->SetPoint(i, x[i], fftinverse->GetPointReal(j)*(fXmax-fXmin)/(fNofPoints*fNofPoints) );
310 }
311 fGraphConv->SetBit(TGraph::kIsSortedX); // indicate that points are sorted in X to speed up TGraph::Eval
312 fFlagGraph = true; // we can use the graph
313
314 // delete the fft objects
315 delete fft1;
316 delete fft2;
317 delete fftinverse;
318}
319
320////////////////////////////////////////////////////////////////////////////////
321
323{
324 if (!fFlagGraph) MakeFFTConv();
325 // if cannot make FFT use numconv
326 if (fGraphConv)
327 return fGraphConv -> Eval(t);
328 else
329
330 return EvalNumConv(t);
331}
332
333////////////////////////////////////////////////////////////////////////////////
334/// Perform numerical convolution.
335/// Could in principle cache the integral in a Graph as it is done for the FFTW
336
338{
339 TF1Convolution_EvalWrapper fconv( *fFunction1, *fFunction2, t);
340 Double_t result = 0;
341
343 if (fXmin != - TMath::Infinity() && fXmax != TMath::Infinity() )
344 result = integrator.Integral(fXmin, fXmax);
345 else if (fXmin == - TMath::Infinity() && fXmax != TMath::Infinity() )
346 result = integrator.IntegralLow(fXmax);
347 else if (fXmin != - TMath::Infinity() && fXmax == TMath::Infinity() )
348 result = integrator.IntegralUp(fXmin);
349 else if (fXmin == - TMath::Infinity() && fXmax == TMath::Infinity() )
350 result = integrator.Integral();
351
352 return result;
353}
354
355////////////////////////////////////////////////////////////////////////////////
356/// Used in TF1 when doing the fit, will be evaluated at each point.
357
359{
360 if (p!=0) TF1Convolution::SetParameters(p); // first refresh the parameters
361
362 Double_t result = 0.;
363 if (fFlagFFT)
364 result = EvalFFTConv(x[0]);
365 else
366 result = EvalNumConv(x[0]);
367 return result;
368}
369
370////////////////////////////////////////////////////////////////////////////////
371
373{
374 if (n<0) return;
375 fNofPoints = n;
376 if (fGraphConv) fGraphConv -> Set(fNofPoints);
377 fFlagGraph = false; // to indicate we need to re-do the graph
378}
379
380////////////////////////////////////////////////////////////////////////////////
381
383{
384 bool equalParams = true;
385 for (int i=0; i<fNofParams1; i++) {
386 fFunction1->SetParameter(i, params[i]);
387 equalParams &= (fParams1[i] == params[i]);
388 fParams1[i] = params[i];
389 }
390 Int_t k = 0;
391 Int_t offset = 0;
392 Int_t offset2 = 0;
393 if (fCstIndex!=-1) offset = 1;
394 Int_t totalnofparams = fNofParams1+fNofParams2+offset;
395 for (int i=fNofParams1; i<totalnofparams; i++) {
396 if (k==fCstIndex)
397 {
398 k++;
399 offset2=1;
400 continue;
401 }
402 fFunction2->SetParameter(k, params[i - offset2]);
403 equalParams &= (fParams2[k - offset2] == params[i - offset2]);
404 fParams2[k - offset2] = params[i - offset2];
405 k++;
406 }
407
408 if (!equalParams) fFlagGraph = false; // to indicate we need to re-do the convolution
409}
410
411////////////////////////////////////////////////////////////////////////////////
412
414 Double_t p4, Double_t p5, Double_t p6, Double_t p7)
415{
416 Double_t params[]={p0,p1,p2,p3,p4,p5,p6,p7};
418}
419
420////////////////////////////////////////////////////////////////////////////////
421
423{
424 if (percentage<0) return;
425 double range = fXmax = fXmin;
426 fXmin -= percentage * range;
427 fXmax += percentage * range;
428 fFlagGraph = false; // to indicate we need to re-do the convolution
429}
430
431////////////////////////////////////////////////////////////////////////////////
432
434{
435 if (a >= b) {
436 Warning("SetRange", "Invalid range: %f >= %f", a, b);
437 return;
438 }
439
440 fXmin = a;
441 fXmax = b;
442 if (fFlagFFT && ( a==-TMath::Infinity() || b==TMath::Infinity() ) )
443 {
444 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.");
445 if (a ==-TMath::Infinity()) fXmin = fFunction1 -> GetXmin();
446 if ( b== TMath::Infinity()) fXmax = fFunction1 -> GetXmax();
447 // add a spill over of 10% in this case
448 SetExtraRange(0.1);
449 }
450 fFlagGraph = false; // to indicate we need to re-do the convolution
451}
452
453////////////////////////////////////////////////////////////////////////////////
454
456{
457 a = fXmin;
458 b = fXmax;
459}
460
461////////////////////////////////////////////////////////////////////////////////
462/// Update the two component functions of the convolution
463
465{
466 fFunction1->Update();
467 fFunction2->Update();
468}
469
470////////////////////////////////////////////////////////////////////////////////
471
473{
474 // copy numbers
475 ((TF1Convolution &)obj).fXmin = fXmin;
476 ((TF1Convolution &)obj).fXmax = fXmax;
477 ((TF1Convolution &)obj).fNofParams1 = fNofParams1;
478 ((TF1Convolution &)obj).fNofParams2 = fNofParams2;
479 ((TF1Convolution &)obj).fCstIndex = fCstIndex;
480 ((TF1Convolution &)obj).fNofPoints = fNofPoints;
481 ((TF1Convolution &)obj).fFlagFFT = fFlagFFT;
482 ((TF1Convolution &)obj).fFlagGraph = false; // since we're not copying the graph
483
484 // copy vectors
485 ((TF1Convolution &)obj).fParams1 = fParams1;
486 ((TF1Convolution &)obj).fParams2 = fParams2;
487 ((TF1Convolution &)obj).fParNames = fParNames;
488
489 // we need to copy the content of the unique_ptr's
490 ((TF1Convolution &)obj).fFunction1 = std::unique_ptr<TF1>((TF1 *)new TF1() );
491 ((TF1Convolution &)obj).fFunction2 = std::unique_ptr<TF1>((TF1 *)new TF1() );
492 fFunction1->Copy(*(((TF1Convolution &)obj).fFunction1 ) );
493 fFunction2->Copy(*(((TF1Convolution &)obj).fFunction2 ) );
494 // fGraphConv is transient anyway, so we don't bother to copy it
495}
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define e(i)
Definition: RSha256.hxx:103
static const double x2[5]
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define ClassImp(name)
Definition: Rtypes.h:365
R__EXTERN Int_t gDebug
Definition: Rtypes.h:91
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
TRObject operator()(const T1 &t1) const
#define gROOT
Definition: TROOT.h:415
static IntegrationOneDim::Type DefaultIntegratorType()
User Class for performing numerical integration of a function in one dimension.
Definition: Integrator.h:94
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
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
double IntegralLow(const IGenFunction &f, double b)
evaluate the Integral of a function f over the over the semi-infinite interval (-inf,...
Definition: Integrator.h:292
Class wrapping convolution of two functions.
std::vector< Double_t > fParams1
void Copy(TObject &obj) const
Copy this to obj.
Int_t fCstIndex
Index of the constant parameter f the first function.
void SetRange(Double_t a, Double_t b)
void GetRange(Double_t &a, Double_t &b) const
std::vector< TString > fParNames
Parameters' names.
std::unique_ptr< TF1 > fFunction1
First function to be convolved.
Int_t fNofPoints
Number of point for FFT array.
Double_t GetXmin() const
std::vector< Double_t > fParams2
Double_t fXmin
Minimal bound of the range of the convolution.
void SetExtraRange(Double_t percentage)
Double_t EvalNumConv(Double_t t)
Perform numerical convolution.
Bool_t fFlagFFT
Choose FFT or numerical convolution.
Double_t GetXmax() const
TF1Convolution & operator=(const TF1Convolution &rhs)
Operator =.
void MakeFFTConv()
Perform the FFT of the two functions.
Double_t operator()(const Double_t *x, const Double_t *p)
Used in TF1 when doing the fit, will be evaluated at each point.
void SetParameters(const Double_t *params)
void SetNofPointsFFT(Int_t n)
TF1Convolution()
constructor without arguments
const char * GetParName(Int_t ipar) const
Double_t EvalFFTConv(Double_t t)
Double_t fXmax
Maximal bound of the range of the convolution.
std::unique_ptr< TF1 > fFunction2
Second function to be convolved.
Bool_t fFlagGraph
! Tells if the graph is already done or not
void InitializeDataMembers(TF1 *function1, TF1 *function2, Bool_t useFFT)
void Update()
Update the two component functions of the convolution.
std::unique_ptr< TGraph > fGraphConv
! Graph of the convolution
1-Dim function class
Definition: TF1.h:211
virtual void Copy(TObject &f1) const
Copy this F1 to a new F1.
Definition: TF1.cxx:995
static void InitStandardFunctions()
Create the basic function objects.
Definition: TF1.cxx:2469
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=0)
Evaluate function with given coordinates and parameters.
Definition: TF1.cxx:1458
@ kNotGlobal
Definition: TF1.h:320
virtual Int_t GetNdim() const
Definition: TF1.h:479
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
@ kIsSortedX
graph is sorted in X points
Definition: TGraph.h:73
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
An array of TObjects.
Definition: TObjArray.h:37
Collectable string class.
Definition: TObjString.h:28
Mother of all ROOT objects.
Definition: TObject.h:37
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
virtual TObject * FindObject(const char *name) const
Must be redefined in derived classes.
Definition: TObject.cxx:321
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:908
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:854
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
TObjArray * Tokenize(const TString &delim) const
This function is used to isolate sequential tokens in a TString.
Definition: TString.cxx:2197
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:2311
TVirtualFFT is an interface class for Fast Fourier Transforms.
Definition: TVirtualFFT.h:88
static TVirtualFFT * FFT(Int_t ndim, Int_t *n, Option_t *option)
Returns a pointer to the FFT of requested size and type.
virtual Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const =0
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
TF1 * f1
Definition: legend1.C:11
Double_t Infinity()
Returns an infinity as defined by the IEEE standard.
Definition: TMath.h:904
auto * a
Definition: textangle.C:12