Logo ROOT   6.10/09
Reference Guide
TF1Convolution.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Authors: L. Moneta, A. Flandi 08/2014
3 //
4 /**********************************************************************
5  * *
6  * Copyright (c) 2015 ROOT Team, CERN/PH-SFT *
7  * *
8  * *
9  **********************************************************************/
10 //
11 // TF1Convolution.cxx
12 //
13 //
14 // Created by AurĂ©lie Flandi on 27.08.14.
15 //
16 //
17 //
18 
19 #include "TF1Convolution.h"
20 #include "Riostream.h"
21 #include "TROOT.h"
22 #include "TObject.h"
23 #include "TObjString.h"
24 #include "TMath.h"
25 #include "Math/Integrator.h"
27 #include "Math/IntegratorOptions.h"
28 #include "Math/GaussIntegrator.h"
31 #include "Math/Functor.h"
32 #include "TVirtualFFT.h"
33 #include "TClass.h"
34 
35 ////////////////////////////////////////////////////////////////////////////////
36 /** \class TF1Convolution
37  \ingroup Hist
38  \brief Class wrapping convolution of two functions
39 
40 Class wrapping convolution of two functions: evaluation of \f$\int f(x)g(x-t)dx\f$
41 
42 The convolution is performed by default using FFTW if it is available .
43 One can pass optionally the range of the convolution (by default the first function range is used).
44 Note that when using Discrete Fouriere Transform (as FFTW), it is a circular transform, so the functions should be
45 approximatly zero at the end of the range. If they are significantly different than zero on one side (e.g. the left side)
46 a spill over will occur on the other side (e.g right side).
47 If no function range is given by default the function1 range + 10% is used
48 One shoud use also a not too small number of points for the DFT (a minimum of 1000). By default 10000 points are used.
49 */////////////////////////////////////////////////////////////////////////////////
50 
51 class TF1Convolution_EvalWrapper
52 {
53  std::shared_ptr < TF1 > fFunction1;
54  std::shared_ptr < TF1 > fFunction2;
55  Double_t fT0;
56 
57 public:
58 
59  TF1Convolution_EvalWrapper(std::shared_ptr<TF1> & f1 , std::shared_ptr<TF1> & f2, Double_t t)
60  : fFunction1(f1), fFunction2(f2), fT0(t)
61  {
62  }
64  {
65  // use EvalPar that is faster
66  Double_t xx[2];
67  xx[0] = x;
68  xx[1] = fT0-x;
69  return fFunction1->EvalPar(xx,nullptr) * fFunction2->EvalPar(xx+1,nullptr);
70  }
71 };
72 
73 ////////////////////////////////////////////////////////////////////////////////
74 /// use copy instead of Clone
75 
76 void TF1Convolution::InitializeDataMembers(TF1* function1, TF1* function2, Bool_t useFFT)
77 {
78  if (function1) {
79  TF1 * fnew1 = (TF1*) function1->IsA()->New();
80  function1->Copy(*fnew1);
81  fFunction1 = std::shared_ptr<TF1>(fnew1);
82  }
83  if (function2) {
84  TF1 * fnew2 = (TF1*) function2->IsA()->New();
85  function2->Copy(*fnew2);
86  fFunction2 = std::shared_ptr<TF1>(fnew2);
87  }
88  if (fFunction1.get() == nullptr|| fFunction2.get() == nullptr)
89  Fatal("InitializeDataMembers","Invalid functions - Abort");
90 
91  // add by default an extra 10% on each side
92  fFunction1->GetRange(fXmin, fXmax);
93  Double_t range = fXmax - fXmin;
94  fXmin -= 0.1*range;
95  fXmax += 0.1*range;
96  fNofParams1 = fFunction1->GetNpar();
97  fNofParams2 = fFunction2->GetNpar();
98  fParams1 = std::vector<Double_t>(fNofParams1);
99  fParams2 = std::vector<Double_t>(fNofParams2);
100  fCstIndex = fFunction2-> GetParNumber("Constant");
101  fFlagFFT = useFFT;
102  fFlagGraph = false;
103  fNofPoints = 10000;
104 
105  //std::cout<<"before: NofParams2 = "<<fNofParams2<<std::endl;
106 
107  fParNames.reserve( fNofParams1 + fNofParams2);
108  for (int i=0; i<fNofParams1; i++)
109  {
110  fParams1[i] = fFunction1 -> GetParameter(i);
111  fParNames.push_back(fFunction1 -> GetParName(i) );
112  }
113  for (int i=0; i<fNofParams2; i++)
114  {
115  fParams2[i] = fFunction2 -> GetParameter(i);
116  if (i != fCstIndex) fParNames.push_back(fFunction2 -> GetParName(i) );
117  }
118  if (fCstIndex!=-1)
119  {
120  fFunction2 -> FixParameter(fCstIndex,1.);
121  fNofParams2 = fNofParams2-1;
122  fParams2.erase(fParams2.begin()+fCstIndex);
123  }
124 }
125 ////////////////////////////////////////////////////////////////////////////////
126 /// constructor from the two function pointer and a flag is using FFT
127 
128 TF1Convolution::TF1Convolution(TF1* function1, TF1* function2, Bool_t useFFT)
129 {
130  InitializeDataMembers(function1,function2, useFFT);
131 }
132 
133 ////////////////////////////////////////////////////////////////////////////////
134 /// constructor from the two function pointer and the convolution range
135 
137 {
138  InitializeDataMembers(function1, function2,useFFT);
139  if (xmin < xmax) {
140  fXmin = xmin;
141  fXmax = xmax;
142  }
143 }
144 
145 ////////////////////////////////////////////////////////////////////////////////
146 /// constructor from a formula expression as f1 * f2 where f1 and f2 are two functions known to ROOT
147 
149 {
151 
152  TObjArray *objarray = formula.Tokenize("*");
153  std::vector < TString > stringarray(2);
154  std::vector < TF1* > funcarray(2);
155  for (int i=0; i<2; i++)
156  {
157  stringarray[i] = ((TObjString*)((*objarray)[i])) -> GetString();
158  stringarray[i].ReplaceAll(" ","");
159  funcarray[i] = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(stringarray[i]));
160  // case function is not found try to use as a TFormula
161  if (funcarray[i] == nullptr) {
162  TF1 * f = new TF1(TString::Format("f_conv_%d",i+1),stringarray[i]);
163  if (!f->GetFormula()->IsValid() )
164  Error("TF1Convolution","Invalid formula : %s",stringarray[i].Data() );
165  if (i == 0)
166  fFunction1 = std::shared_ptr<TF1>(f);
167  else
168  fFunction2 = std::shared_ptr<TF1>(f);
169  }
170  }
171  InitializeDataMembers(funcarray[0], funcarray[1],useFFT);
172  if (xmin < xmax) {
173  fXmin = xmin;
174  fXmax = xmax;
175  }
176 }
177 
178 ////////////////////////////////////////////////////////////////////////////////
179 /// constructor from 2 function names where f1 and f2 are two functions known to ROOT
180 /// if the function names are not knwon to ROOT then a corresponding
181 
183 {
185  (TString)formula1.ReplaceAll(" ","");
186  (TString)formula2.ReplaceAll(" ","");
187  TF1* f1 = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(formula1));
188  TF1* f2 = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(formula2));
189  // if function do not exists try using TFormula
190  if (!f1) {
191  fFunction1 = std::shared_ptr<TF1>(new TF1("f_conv_1",formula1) );
192  if (!fFunction1->GetFormula()->IsValid() )
193  Error("TF1Convolution","Invalid formula for : %s",formula1.Data() );
194  }
195  if (!f2) {
196  fFunction2 = std::shared_ptr<TF1>(new TF1("f_conv_1",formula2) );
197  if (!fFunction2->GetFormula()->IsValid() )
198  Error("TF1Convolution","Invalid formula for : %s",formula2.Data() );
199  }
200  // if f1 or f2 are null ptr are not used in InitializeDataMembers
201  InitializeDataMembers(f1, f2,useFFT);
202  if (xmin < xmax) {
203  fXmin = xmin;
204  fXmax = xmax;
205  }
206 }
207 
208 ////////////////////////////////////////////////////////////////////////////////
209 ///perform the FFT of the two functions
210 
212 {
213  if (gDebug)
214  Info("MakeFFTConv","Making FFT convolution using %d points in range [%g,%g]",fNofPoints,fXmin,fXmax);
215 
216  std::vector < Double_t > x (fNofPoints);
217  std::vector < Double_t > in1(fNofPoints);
218  std::vector < Double_t > in2(fNofPoints);
219 
220  TVirtualFFT *fft1 = TVirtualFFT::FFT(1, &fNofPoints, "R2C K");
221  TVirtualFFT *fft2 = TVirtualFFT::FFT(1, &fNofPoints, "R2C K");
222  if (fft1 == nullptr || fft2 == nullptr) {
223  Warning("MakeFFTConv","Cannot use FFT, probably FFTW package is not available. Switch to numerical convolution");
224  fFlagFFT = false;
225  return;
226  }
227 
228  // apply a shift in order to have the second function centered around middle of the range of the convolution
229  Double_t shift2 = 0.5*(fXmin+fXmax);
230  Double_t x2;
231  for (int i=0; i<fNofPoints; i++)
232  {
233  x[i] = fXmin + (fXmax-fXmin)/(fNofPoints-1)*i;
234  x2 = x[i] - shift2;
235  in1[i] = fFunction1 -> EvalPar( &x[i], nullptr);
236  in2[i] = fFunction2 -> EvalPar( &x2, nullptr);
237  fft1 -> SetPoint(i, in1[i]);
238  fft2 -> SetPoint(i, in2[i]);
239  }
240  fft1 -> Transform();
241  fft2 -> Transform();
242 
243  //inverse transformation of the product
244 
245  TVirtualFFT *fftinverse = TVirtualFFT::FFT(1, &fNofPoints, "C2R K");
246  Double_t re1, re2, im1, im2, out_re, out_im;
247 
248  for (int i=0;i<=fNofPoints/2.;i++)
249  {
250  fft1 -> GetPointComplex(i,re1,im1);
251  fft2 -> GetPointComplex(i,re2,im2);
252  out_re = re1*re2 - im1*im2;
253  out_im = re1*im2 + re2*im1;
254  fftinverse -> SetPoint(i, out_re, out_im);
255  }
256  fftinverse -> Transform();
257 
258  // fill a graph with the result of the convolution
259  if (!fGraphConv) fGraphConv = std::shared_ptr< TGraph >(new TGraph(fNofPoints));
260 
261  for (int i=0;i<fNofPoints;i++)
262  {
263  // we need this since we have applied a shift in the middle of f2
264  int j = i + fNofPoints/2;
265  if (j >= fNofPoints) j -= fNofPoints;
266  // need to normalize by dividing by the number of points and multiply by the bin width = Range/Number of points
267  fGraphConv->SetPoint(i, x[i], fftinverse->GetPointReal(j)*(fXmax-fXmin)/(fNofPoints*fNofPoints) );
268  }
269  fGraphConv->SetBit(TGraph::kIsSortedX); // indicate that points are sorted in X to speed up TGraph::Eval
270  fFlagGraph = true; // we can use the graph
271 }
272 
273 ////////////////////////////////////////////////////////////////////////////////
274 
276 {
277  if (!fFlagGraph) MakeFFTConv();
278  // if cannot make FFT use numconv
279  if (fGraphConv)
280  return fGraphConv -> Eval(t);
281  else
282  return EvalNumConv(t);
283 }
284 
285 ////////////////////////////////////////////////////////////////////////////////
286 /// perform numerical convolution
287 /// could in principle cache the integral in a Graph as it is done for the FFTW
288 
290 {
291  TF1Convolution_EvalWrapper fconv( fFunction1, fFunction2, t);
292  Double_t result = 0;
293 
295  if (fXmin != - TMath::Infinity() && fXmax != TMath::Infinity() )
296  result = integrator.Integral(fXmin, fXmax);
297  else if (fXmin == - TMath::Infinity() && fXmax != TMath::Infinity() )
298  result = integrator.IntegralLow(fXmax);
299  else if (fXmin != - TMath::Infinity() && fXmax == TMath::Infinity() )
300  result = integrator.IntegralUp(fXmin);
301  else if (fXmin == - TMath::Infinity() && fXmax == TMath::Infinity() )
302  result = integrator.Integral();
303 
304  return result;
305 }
306 
307 ////////////////////////////////////////////////////////////////////////////////
308 
309 Double_t TF1Convolution::operator()(Double_t* t, Double_t* p)//used in TF1 when doing the fit, will be valuated at each point
310 {
311  if (p!=0) TF1Convolution::SetParameters(p); // first refresh the parameters
312 
313  Double_t result = 0.;
314  if (fFlagFFT) result = EvalFFTConv(t[0]);
315  else result = EvalNumConv(t[0]);
316  return result;
317 }
318 ////////////////////////////////////////////////////////////////////////////////
319 
321 {
322  if (n<0) return;
323  fNofPoints = n;
324  if (fGraphConv) fGraphConv -> Set(fNofPoints); //set nof points of the Tgraph
325  fFlagGraph = false; // to indicate we need to re-do the graph
326 }
327 
328 ////////////////////////////////////////////////////////////////////////////////
329 
331 {
332  bool equalParams = true;
333  for (int i=0; i<fNofParams1; i++) {
334  fFunction1 -> SetParameter(i,p[i]);
335  equalParams &= ( fParams1[i] == p[i] );
336  fParams1[i] = p[i];
337  }
338  Int_t k = 0;
339  Int_t offset = 0;
340  Int_t offset2 = 0;
341  if (fCstIndex!=-1) offset = 1;
342  Int_t totalnofparams = fNofParams1+fNofParams2+offset;
343  for (int i=fNofParams1; i<totalnofparams; i++) {
344  if (k==fCstIndex)
345  {
346  k++;
347  offset2=1;
348  continue;
349  }
350  fFunction2 -> SetParameter(k,p[i-offset2]);
351  equalParams &= ( fParams2[k-offset2] == p[i-offset2] );
352  fParams2[k-offset2] = p[i-offset2];
353  k++;
354  }
355  // std::cout << "parameters for function1 : ";
356  // for (int i = 0; i < fFunction1->GetNpar(); ++i)
357  // std::cout << fFunction1->GetParameter(i) << " ";
358  // std::cout << "\nparameters for function2 : ";
359  // for (int i = 0; i < fFunction2->GetNpar(); ++i)
360  // std::cout << fFunction2->GetParameter(i) << " ";
361  // std::cout << std::endl;
362 
363 //do the graph for FFT convolution
364  if (!equalParams) fFlagGraph = false; // to indicate we need to re-do the convolution
365  // if (fFlagFFT)
366  // {
367  // MakeFFTConv();
368  // }
369 }
370 
371 ////////////////////////////////////////////////////////////////////////////////
372 
374  Double_t p4, Double_t p5, Double_t p6, Double_t p7)
375 {
376  Double_t params[]={p0,p1,p2,p3,p4,p5,p6,p7};
378 }
379 
380 ////////////////////////////////////////////////////////////////////////////////
381 
383 {
384  if (percentage<0) return;
385  double range = fXmax = fXmin;
386  fXmin -= percentage * range;
387  fXmax += percentage * range;
388  fFlagGraph = false; // to indicate we need to re-do the convolution
389 }
390 
391 ////////////////////////////////////////////////////////////////////////////////
392 
394 {
395  if (a>=b) return;
396  fXmin = a;
397  fXmax = b;
398  if (fFlagFFT && ( a==-TMath::Infinity() || b==TMath::Infinity() ) )
399  {
400  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.");
401  if (a ==-TMath::Infinity()) fXmin = fFunction1 -> GetXmin();
402  if ( b== TMath::Infinity()) fXmax = fFunction1 -> GetXmax();
403  // add a spill over of 10% in this case
404  SetExtraRange(0.1);
405  }
406  fFlagGraph = false; // to indicate we need to re-do the convolution
407 }
An array of TObjects.
Definition: TObjArray.h:37
float xmin
Definition: THbookFile.cxx:93
virtual TFormula * GetFormula()
Definition: TF1.h:407
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:775
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:640
#define gROOT
Definition: TROOT.h:375
Basic string class.
Definition: TString.h:129
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
TRObject operator()(const T1 &t1) const
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:2345
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,...)
std::vector< std::vector< double > > Data
void InitializeDataMembers(TF1 *function1, TF1 *function2, Bool_t useFFT)
use copy instead of Clone
static IntegrationOneDim::Type DefaultIntegratorType()
Double_t Infinity()
Definition: TMath.h:797
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
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 could in principle cache the integral in a Graph as it is done for the ...
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:2276
static double p1(double t, double a, double b)
float xmax
Definition: THbookFile.cxx:93
TF1Convolution(TF1 *function1, TF1 *function2, Bool_t useFFT=true)
constructor from the two function pointer and a flag is using FFT
void Warning(const char *location, const char *msgfmt,...)
graph is sorted in X points
Definition: TGraph.h:73
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:2251
double f(double x)
double Double_t
Definition: RtypesCore.h:55
Double_t operator()(Double_t *t, Double_t *p)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
void SetParameters(Double_t *p)
Double_t EvalFFTConv(Double_t t)
double f2(const double *x)
1-Dim function class
Definition: TF1.h:150
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:83
double result[121]
Bool_t IsValid() const
Definition: TFormula.h:181
const Int_t n
Definition: legend1.C:16
void SetExtraRange(Double_t percentage)
const char * Data() const
Definition: TString.h:347