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