Logo ROOT  
Reference Guide
TF3.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Author: Rene Brun 27/10/95
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #include "TROOT.h"
13 #include "TF3.h"
14 #include "TBuffer.h"
15 #include "TMath.h"
16 #include "TH3.h"
17 #include "TVirtualPad.h"
18 #include "TRandom.h"
19 #include "TVectorD.h"
20 #include "Riostream.h"
21 #include "TColor.h"
22 #include "TVirtualFitter.h"
23 #include "TVirtualHistPainter.h"
24 #include "Math/IntegratorOptions.h"
25 #include <cassert>
26 
27 ClassImp(TF3);
28 
29 /** \class TF3
30  \ingroup Hist
31 A 3-Dim function with parameters
32 */
33 
34 ////////////////////////////////////////////////////////////////////////////////
35 /// F3 default constructor
36 
38 {
39  fNpz = 0;
40  fZmin = 0;
41  fZmax = 1;
42 }
43 
44 
45 ////////////////////////////////////////////////////////////////////////////////
46 /// F3 constructor using a formula definition
47 ///
48 /// See TFormula constructor for explanation of the formula syntax.
49 
50 TF3::TF3(const char *name,const char *formula, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Double_t zmin, Double_t zmax, Option_t * opt)
51  :TF2(name,formula,xmin,xmax,ymax,ymin,opt)
52 {
53  fZmin = zmin;
54  fZmax = zmax;
55  fNpz = 30;
56  Int_t ndim = GetNdim();
57  // accept 1-d or 2-d formula
58  if (ndim < 3) fNdim = 3;
59  if (ndim > 3 && xmin < xmax && ymin < ymax && zmin < zmax) {
60  Error("TF3","function: %s/%s has dimension %d instead of 3",name,formula,ndim);
61  MakeZombie();
62  }
63 }
64 
65 ////////////////////////////////////////////////////////////////////////////////
66 /// F3 constructor using a pointer to real function
67 ///
68 /// \param[in] npar is the number of free parameters used by the function
69 ///
70 /// For example, for a 3-dim function with 3 parameters, the user function
71 /// looks like:
72 ///
73 /// Double_t fun1(Double_t *x, Double_t *par)
74 /// return par[0]*x[2] + par[1]*exp(par[2]*x[0]*x[1]);
75 ///
76 /// WARNING! A function created with this constructor cannot be Cloned.
77 
79  :TF2(name,fcn,xmin,xmax,ymin,ymax,npar,ndim)
80 {
81  fZmin = zmin;
82  fZmax = zmax;
83  fNpz = 30;
84 }
85 
86 ////////////////////////////////////////////////////////////////////////////////
87 /// F3 constructor using a pointer to real function---
88 ///
89 /// \param[in] npar is the number of free parameters used by the function
90 ///
91 /// For example, for a 3-dim function with 3 parameters, the user function
92 /// looks like:
93 ///
94 /// Double_t fun1(Double_t *x, Double_t *par)
95 /// return par[0]*x[2] + par[1]*exp(par[2]*x[0]*x[1]);
96 ///
97 /// WARNING! A function created with this constructor cannot be Cloned.
98 
99 TF3::TF3(const char *name,Double_t (*fcn)(const Double_t *, const Double_t *), Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Double_t zmin, Double_t zmax, Int_t npar, Int_t ndim)
100  : TF2(name,fcn,xmin,xmax,ymin,ymax,npar,ndim),
101  fZmin(zmin),
102  fZmax(zmax),
103  fNpz(30)
104 {
105 }
106 
107 ////////////////////////////////////////////////////////////////////////////////
108 /// F3 constructor using a ParamFunctor
109 ///
110 /// a functor class implementing operator() (double *, double *)
111 ///
112 /// \param[in] npar is the number of free parameters used by the function
113 ///
114 /// WARNING! A function created with this constructor cannot be Cloned.
115 
117  : TF2(name, f, xmin, xmax, ymin, ymax, npar, ndim),
118  fZmin(zmin),
119  fZmax(zmax),
120  fNpz(30)
121 {
122 }
123 
124 ////////////////////////////////////////////////////////////////////////////////
125 /// Operator =
126 
127 TF3& TF3::operator=(const TF3 &rhs)
128 {
129  if (this != &rhs) {
130  rhs.Copy(*this);
131  }
132  return *this;
133 }
134 
135 ////////////////////////////////////////////////////////////////////////////////
136 /// F3 default destructor
137 
139 {
140 }
141 
142 ////////////////////////////////////////////////////////////////////////////////
143 /// Copy constructor.
144 
145 TF3::TF3(const TF3 &f3) : TF2()
146 {
147  ((TF3&)f3).Copy(*this);
148 }
149 
150 ////////////////////////////////////////////////////////////////////////////////
151 /// Copy this F3 to a new F3
152 
153 void TF3::Copy(TObject &obj) const
154 {
155  TF2::Copy(obj);
156  ((TF3&)obj).fZmin = fZmin;
157  ((TF3&)obj).fZmax = fZmax;
158  ((TF3&)obj).fNpz = fNpz;
159 }
160 
161 ////////////////////////////////////////////////////////////////////////////////
162 /// Compute distance from point px,py to a function
163 ///
164 /// Compute the closest distance of approach from point px,py to this function.
165 /// The distance is computed in pixels units.
166 
167 
169 {
170  return TF1::DistancetoPrimitive(px, py);
171 
172 }
173 
174 ////////////////////////////////////////////////////////////////////////////////
175 /// Draw this function with its current attributes
176 
177 void TF3::Draw(Option_t *option)
178 {
179  TString opt = option;
180  opt.ToLower();
181  if (gPad && !opt.Contains("same")) gPad->Clear();
182 
183  AppendPad(option);
184 
185 }
186 
187 ////////////////////////////////////////////////////////////////////////////////
188 /// Execute action corresponding to one event
189 ///
190 /// This member function is called when a F3 is clicked with the locator
191 
192 void TF3::ExecuteEvent(Int_t event, Int_t px, Int_t py)
193 {
194  TF1::ExecuteEvent(event, px, py);
195 }
196 
197 ////////////////////////////////////////////////////////////////////////////////
198 /// Return minimum/maximum value of the function
199 ///
200 /// To find the minimum on a range, first set this range via the SetRange function
201 /// If a vector x of coordinate is passed it will be used as starting point for the minimum.
202 /// In addition on exit x will contain the coordinate values at the minimuma
203 /// If x is NULL or x is inifinity or NaN, first, a grid search is performed to find the initial estimate of the
204 /// minimum location. The range of the function is divided into fNpx and fNpy
205 /// sub-ranges. If the function is "good" (or "bad"), these values can be changed
206 /// by SetNpx and SetNpy functions
207 ///
208 /// Then, a minimization is used with starting values found by the grid search
209 /// The minimizer algorithm used (by default Minuit) can be changed by callinga
210 /// ROOT::Math::Minimizer::SetDefaultMinimizerType("..")
211 /// Other option for the minimizer can be set using the static method of the MinimizerOptions class
212 
214 {
215  //First do a grid search with step size fNpx and fNpy
216 
217  Double_t xx[3];
218  Double_t rsign = (findmax) ? -1. : 1.;
219  TF3 & function = const_cast<TF3&>(*this); // needed since EvalPar is not const
220  Double_t xxmin = 0, yymin = 0, zzmin = 0, ttmin = 0;
221  if (x == NULL || ( (x!= NULL) && ( !TMath::Finite(x[0]) || !TMath::Finite(x[1]) || !TMath::Finite(x[2]) ) ) ){
222  Double_t dx = (fXmax - fXmin)/fNpx;
223  Double_t dy = (fYmax - fYmin)/fNpy;
224  Double_t dz = (fZmax - fZmin)/fNpz;
225  xxmin = fXmin;
226  yymin = fYmin;
227  zzmin = fZmin;
228  ttmin = rsign * TMath::Infinity();
229  for (Int_t i=0; i<fNpx; i++){
230  xx[0]=fXmin + (i+0.5)*dx;
231  for (Int_t j=0; j<fNpy; j++){
232  xx[1]=fYmin+(j+0.5)*dy;
233  for (Int_t k=0; k<fNpz; k++){
234  xx[2] = fZmin+(k+0.5)*dz;
235  Double_t tt = function(xx);
236  if (rsign*tt < rsign*ttmin) {xxmin = xx[0], yymin = xx[1]; zzmin = xx[2]; ttmin=tt;}
237  }
238  }
239  }
240 
241  xxmin = TMath::Min(fXmax, xxmin);
242  yymin = TMath::Min(fYmax, yymin);
243  zzmin = TMath::Min(fZmax, zzmin);
244  }
245  else {
246  xxmin = x[0];
247  yymin = x[1];
248  zzmin = x[2];
249  zzmin = function(xx);
250  }
251  xx[0] = xxmin;
252  xx[1] = yymin;
253  xx[2] = zzmin;
254 
255  double fmin = GetMinMaxNDim(xx,findmax);
256  if (rsign*fmin < rsign*zzmin) {
257  if (x) {x[0] = xx[0]; x[1] = xx[1]; x[2] = xx[2];}
258  return fmin;
259  }
260  // here if minimization failed
261  if (x) { x[0] = xxmin; x[1] = yymin; x[2] = zzmin; }
262  return ttmin;
263 }
264 
265 ////////////////////////////////////////////////////////////////////////////////
266 /// Compute the X, Y and Z values corresponding to the minimum value of the function
267 /// on its range.
268 ///
269 /// Returns the function value at the minimum.
270 /// To find the minimum on a subrange, use the SetRange() function first.
271 ///
272 /// Method:
273 /// First, a grid search is performed to find the initial estimate of the
274 /// minimum location. The range of the function is divided
275 /// into fNpx,fNpy and fNpz sub-ranges. If the function is "good" (or "bad"),
276 /// these values can be changed by SetNpx(), SetNpy() and SetNpz() functions.
277 /// Then, Minuit minimization is used with starting values found by the grid search
278 ///
279 /// Note that this method will always do first a grid search in contrast to GetMinimum
280 
282 {
283  double xx[3] = { 0,0,0 };
284  xx[0] = TMath::QuietNaN(); // to force to do grid search in TF3::FindMinMax
285  double fmin = FindMinMax(xx, false);
286  x = xx[0]; y = xx[1]; z = xx[2];
287  return fmin;
288 
289 }
290 
291 ////////////////////////////////////////////////////////////////////////////////
292 /// Compute the X, Y and Z values corresponding to the maximum value of the function
293 /// on its range.
294 ///
295 /// Return the function value at the maximum. See TF3::GetMinimumXYZ
296 
298 {
299  double xx[3] = { 0,0,0 };
300  xx[0] = TMath::QuietNaN(); // to force to do grid search in TF3::FindMinMax
301  double fmax = FindMinMax(xx, true);
302  x = xx[0]; y = xx[1]; z = xx[2];
303  return fmax;
304 
305 }
306 
307 ////////////////////////////////////////////////////////////////////////////////
308 /// Return 3 random numbers following this function shape
309 ///
310 /// The distribution contained in this TF3 function is integrated
311 /// over the cell contents.
312 /// It is normalized to 1.
313 /// Getting the three random numbers implies:
314 /// - Generating a random number between 0 and 1 (say r1)
315 /// - Look in which cell in the normalized integral r1 corresponds to
316 /// - make a linear interpolation in the returned cell
317 ///
318 /// IMPORTANT NOTE
319 ///
320 /// The integral of the function is computed at fNpx * fNpy * fNpz points.
321 /// If the function has sharp peaks, you should increase the number of
322 /// points (SetNpx, SetNpy, SetNpz) such that the peak is correctly tabulated
323 /// at several points.
324 
325 void TF3::GetRandom3(Double_t &xrandom, Double_t &yrandom, Double_t &zrandom, TRandom * rng)
326 {
327  // Check if integral array must be build
328  Int_t i,j,k,cell;
329  Double_t dx = (fXmax-fXmin)/fNpx;
330  Double_t dy = (fYmax-fYmin)/fNpy;
331  Double_t dz = (fZmax-fZmin)/fNpz;
332  Int_t ncells = fNpx*fNpy*fNpz;
333  Double_t xx[3];
334  Double_t *parameters = GetParameters();
335  InitArgs(xx,parameters);
336  if (fIntegral.empty() ) {
337  fIntegral.resize(ncells+1);
338  //fIntegral = new Double_t[ncells+1];
339  fIntegral[0] = 0;
340  Double_t integ;
341  Int_t intNegative = 0;
342  cell = 0;
343  for (k=0;k<fNpz;k++) {
344  xx[2] = fZmin+(k+0.5)*dz;
345  for (j=0;j<fNpy;j++) {
346  xx[1] = fYmin+(j+0.5)*dy;
347  for (i=0;i<fNpx;i++) {
348  xx[0] = fXmin+(i+0.5)*dx;
349  integ = EvalPar(xx,parameters);
350  if (integ < 0) {intNegative++; integ = -integ;}
351  fIntegral[cell+1] = fIntegral[cell] + integ;
352  cell++;
353  }
354  }
355  }
356  if (intNegative > 0) {
357  Warning("GetRandom3","function:%s has %d negative values: abs assumed",GetName(),intNegative);
358  }
359  if (fIntegral[ncells] == 0) {
360  Error("GetRandom3","Integral of function is zero");
361  return;
362  }
363  for (i=1;i<=ncells;i++) { // normalize integral to 1
364  fIntegral[i] /= fIntegral[ncells];
365  }
366  }
367 
368 // return random numbers
369  Double_t r;
370  if (!rng) rng = gRandom;
371  r = rng->Rndm();
372  cell = TMath::BinarySearch(ncells,fIntegral.data(),r);
373  k = cell/(fNpx*fNpy);
374  j = (cell -k*fNpx*fNpy)/fNpx;
375  i = cell -fNpx*(j +fNpy*k);
376  xrandom = fXmin +dx*i +dx*rng->Rndm();
377  yrandom = fYmin +dy*j +dy*rng->Rndm();
378  zrandom = fZmin +dz*k +dz*rng->Rndm();
379 }
380 
381 ////////////////////////////////////////////////////////////////////////////////
382 /// Return range of function
383 
385 {
386  xmin = fXmin;
387  xmax = fXmax;
388  ymin = fYmin;
389  ymax = fYmax;
390  zmin = fZmin;
391  zmax = fZmax;
392 }
393 
394 
395 ////////////////////////////////////////////////////////////////////////////////
396 /// Get value corresponding to X in array of fSave values
397 
399 {
400  //if (fNsave <= 0) return 0;
401  if (fSave.empty()) return 0;
402  Int_t np = fSave.size() - 9;
403  Double_t xmin = Double_t(fSave[np+0]);
404  Double_t xmax = Double_t(fSave[np+1]);
405  Double_t ymin = Double_t(fSave[np+2]);
406  Double_t ymax = Double_t(fSave[np+3]);
407  Double_t zmin = Double_t(fSave[np+4]);
408  Double_t zmax = Double_t(fSave[np+5]);
409  Int_t npx = Int_t(fSave[np+6]);
410  Int_t npy = Int_t(fSave[np+7]);
411  Int_t npz = Int_t(fSave[np+8]);
412  Double_t x = Double_t(xx[0]);
413  Double_t dx = (xmax-xmin)/npx;
414  if (x < xmin || x > xmax) return 0;
415  if (dx <= 0) return 0;
416  Double_t y = Double_t(xx[1]);
417  Double_t dy = (ymax-ymin)/npy;
418  if (y < ymin || y > ymax) return 0;
419  if (dy <= 0) return 0;
420  Double_t z = Double_t(xx[2]);
421  Double_t dz = (zmax-zmin)/npz;
422  if (z < zmin || z > zmax) return 0;
423  if (dz <= 0) return 0;
424 
425  //we make a trilinear interpolation using the 8 points surrounding x,y,z
426  Int_t ibin = Int_t((x-xmin)/dx);
427  Int_t jbin = Int_t((y-ymin)/dy);
428  Int_t kbin = Int_t((z-zmin)/dz);
429  Double_t xlow = xmin + ibin*dx;
430  Double_t ylow = ymin + jbin*dy;
431  Double_t zlow = zmin + kbin*dz;
432  Double_t t = (x-xlow)/dx;
433  Double_t u = (y-ylow)/dy;
434  Double_t v = (z-zlow)/dz;
435  Int_t k1 = (ibin ) + (npx+1)*((jbin ) + (npy+1)*(kbin ));
436  Int_t k2 = (ibin+1) + (npx+1)*((jbin ) + (npy+1)*(kbin ));
437  Int_t k3 = (ibin+1) + (npx+1)*((jbin+1) + (npy+1)*(kbin ));
438  Int_t k4 = (ibin ) + (npx+1)*((jbin+1) + (npy+1)*(kbin ));
439  Int_t k5 = (ibin ) + (npx+1)*((jbin ) + (npy+1)*(kbin+1));
440  Int_t k6 = (ibin+1) + (npx+1)*((jbin ) + (npy+1)*(kbin+1));
441  Int_t k7 = (ibin+1) + (npx+1)*((jbin+1) + (npy+1)*(kbin+1));
442  Int_t k8 = (ibin ) + (npx+1)*((jbin+1) + (npy+1)*(kbin+1));
443  Double_t r = (1-t)*(1-u)*(1-v)*fSave[k1] + t*(1-u)*(1-v)*fSave[k2] + t*u*(1-v)*fSave[k3] + (1-t)*u*(1-v)*fSave[k4] +
444  (1-t)*(1-u)*v*fSave[k5] + t*(1-u)*v*fSave[k6] + t*u*v*fSave[k7] + (1-t)*u*v*fSave[k8];
445  return r;
446 }
447 
448 ////////////////////////////////////////////////////////////////////////////////
449 /// Return Integral of a 3d function in range [ax,bx],[ay,by],[az,bz]
450 /// with a desired relative accuracy.
451 
453 {
454  Double_t a[3], b[3];
455  a[0] = ax;
456  b[0] = bx;
457  a[1] = ay;
458  b[1] = by;
459  a[2] = az;
460  b[2] = bz;
461  Double_t relerr = 0;
462  Int_t n = 3;
464  Int_t nfnevl, ifail;
465  Double_t result = IntegralMultiple(n,a,b,maxpts,epsrel,epsrel, relerr,nfnevl,ifail);
466  if (ifail > 0) {
467  Warning("Integral","failed for %s code=%d, maxpts=%d, epsrel=%g, nfnevl=%d, relerr=%g ",GetName(),ifail,maxpts,epsrel,nfnevl,relerr);
468  }
469  if (gDebug) {
470  Info("Integral","Integral of %s using %d and tol=%f is %f , relerr=%f nfcn=%d",GetName(),maxpts,epsrel,result,relerr,nfnevl);
471  }
472  return result;
473 }
474 
475 ////////////////////////////////////////////////////////////////////////////////
476 /// Return kTRUE is the point is inside the function range
477 
479 {
480  if (x[0] < fXmin || x[0] > fXmax) return kFALSE;
481  if (x[1] < fYmin || x[1] > fYmax) return kFALSE;
482  if (x[2] < fZmin || x[2] > fZmax) return kFALSE;
483  return kTRUE;
484 }
485 
486 ////////////////////////////////////////////////////////////////////////////////
487 /// Create a histogram for axis range.
488 
490 {
491  TH1* h = new TH3F("R__TF3",(char*)GetTitle(),fNpx,fXmin,fXmax
492  ,fNpy,fYmin,fYmax
493  ,fNpz,fZmin,fZmax);
494  h->SetDirectory(0);
495  return h;
496 }
497 
498 ////////////////////////////////////////////////////////////////////////////////
499 /// Paint this 3-D function with its current attributes
500 
501 void TF3::Paint(Option_t *option)
502 {
503 
504  TString opt = option;
505  opt.ToLower();
506 
507 //- Create a temporary histogram and fill each channel with the function value
508  if (!fHistogram) {
509  fHistogram = new TH3F("R__TF3",(char*)GetTitle(),fNpx,fXmin,fXmax
510  ,fNpy,fYmin,fYmax
511  ,fNpz,fZmin,fZmax);
513  }
514 
515  fHistogram->GetPainter(option)->ProcessMessage("SetF3",this);
516 
517  if (opt.Index("tf3") == kNPOS)
518  opt.Append("tf3");
519 
520  fHistogram->Paint(opt.Data());
521 }
522 
523 ////////////////////////////////////////////////////////////////////////////////
524 /// Set the function clipping box (for drawing) "off".
525 
527 {
528  fClipBoxOn = kFALSE;
529  fClipBox[0] = fClipBox[1] = fClipBox[2] = 0;
530 }
531 
532 ////////////////////////////////////////////////////////////////////////////////
533 /// Save values of function in array fSave
534 
536 {
537  if (!fSave.empty()) fSave.clear();
538  Int_t nsave = (fNpx+1)*(fNpy+1)*(fNpz+1);
539  Int_t fNsave = nsave+9;
540  assert(fNsave > 9);
541  //fSave = new Double_t[fNsave];
542  fSave.resize(fNsave);
543  Int_t i,j,k,l=0;
544  Double_t dx = (xmax-xmin)/fNpx;
545  Double_t dy = (ymax-ymin)/fNpy;
546  Double_t dz = (zmax-zmin)/fNpz;
547  if (dx <= 0) {
548  dx = (fXmax-fXmin)/fNpx;
549  xmin = fXmin +0.5*dx;
550  xmax = fXmax -0.5*dx;
551  }
552  if (dy <= 0) {
553  dy = (fYmax-fYmin)/fNpy;
554  ymin = fYmin +0.5*dy;
555  ymax = fYmax -0.5*dy;
556  }
557  if (dz <= 0) {
558  dz = (fZmax-fZmin)/fNpz;
559  zmin = fZmin +0.5*dz;
560  zmax = fZmax -0.5*dz;
561  }
562  Double_t xv[3];
563  Double_t *pp = GetParameters();
564  InitArgs(xv,pp);
565  for (k=0;k<=fNpz;k++) {
566  xv[2] = zmin + dz*k;
567  for (j=0;j<=fNpy;j++) {
568  xv[1] = ymin + dy*j;
569  for (i=0;i<=fNpx;i++) {
570  xv[0] = xmin + dx*i;
571  fSave[l] = EvalPar(xv,pp);
572  l++;
573  }
574  }
575  }
576  fSave[nsave+0] = xmin;
577  fSave[nsave+1] = xmax;
578  fSave[nsave+2] = ymin;
579  fSave[nsave+3] = ymax;
580  fSave[nsave+4] = zmin;
581  fSave[nsave+5] = zmax;
582  fSave[nsave+6] = fNpx;
583  fSave[nsave+7] = fNpy;
584  fSave[nsave+8] = fNpz;
585 }
586 
587 ////////////////////////////////////////////////////////////////////////////////
588 /// Save primitive as a C++ statement(s) on output stream out
589 
590 void TF3::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
591 {
592  char quote = '"';
593  out<<" "<<std::endl;
594  if (gROOT->ClassSaved(TF3::Class())) {
595  out<<" ";
596  } else {
597  out<<" TF3 *";
598  }
599  if (!fMethodCall) {
600  out<<GetName()<<" = new TF3("<<quote<<GetName()<<quote<<","<<quote<<GetTitle()<<quote<<","<<fXmin<<","<<fXmax<<","<<fYmin<<","<<fYmax<<","<<fZmin<<","<<fZmax<<");"<<std::endl;
601  } else {
602  out<<GetName()<<" = new TF3("<<quote<<GetName()<<quote<<","<<GetTitle()<<","<<fXmin<<","<<fXmax<<","<<fYmin<<","<<fYmax<<","<<fZmin<<","<<fZmax<<","<<GetNpar()<<");"<<std::endl;
603  }
604 
605  if (GetFillColor() != 0) {
606  if (GetFillColor() > 228) {
608  out<<" "<<GetName()<<"->SetFillColor(ci);" << std::endl;
609  } else
610  out<<" "<<GetName()<<"->SetFillColor("<<GetFillColor()<<");"<<std::endl;
611  }
612  if (GetLineColor() != 1) {
613  if (GetLineColor() > 228) {
615  out<<" "<<GetName()<<"->SetLineColor(ci);" << std::endl;
616  } else
617  out<<" "<<GetName()<<"->SetLineColor("<<GetLineColor()<<");"<<std::endl;
618  }
619  if (GetNpz() != 100) {
620  out<<" "<<GetName()<<"->SetNpz("<<GetNpz()<<");"<<std::endl;
621  }
622  if (GetChisquare() != 0) {
623  out<<" "<<GetName()<<"->SetChisquare("<<GetChisquare()<<");"<<std::endl;
624  }
625  Double_t parmin, parmax;
626  for (Int_t i=0;i<GetNpar();i++) {
627  out<<" "<<GetName()<<"->SetParameter("<<i<<","<<GetParameter(i)<<");"<<std::endl;
628  out<<" "<<GetName()<<"->SetParError("<<i<<","<<GetParError(i)<<");"<<std::endl;
629  GetParLimits(i,parmin,parmax);
630  out<<" "<<GetName()<<"->SetParLimits("<<i<<","<<parmin<<","<<parmax<<");"<<std::endl;
631  }
632  out<<" "<<GetName()<<"->Draw("
633  <<quote<<option<<quote<<");"<<std::endl;
634 }
635 
636 ////////////////////////////////////////////////////////////////////////////////
637 /// Set the function clipping box (for drawing) "on" and define the clipping box.
638 /// xclip, yclip and zclip is a point within the function range. All the
639 /// function values having x<=xclip and y<=yclip and z>=zclip are clipped.
640 
642 {
643  fClipBoxOn = kTRUE;
644  fClipBox[0] = xclip;
645  fClipBox[1] = yclip;
646  fClipBox[2] = zclip;
647 }
648 
649 ////////////////////////////////////////////////////////////////////////////////
650 /// Set the number of points used to draw the function
651 ///
652 /// The default number of points along x is 30 for 2-d/3-d functions.
653 /// You can increase this value to get a better resolution when drawing
654 /// pictures with sharp peaks or to get a better result when using TF3::GetRandom2
655 /// the minimum number of points is 4, the maximum is 10000 for 2-d/3-d functions
656 
658 {
659  if (npz < 4) {
660  Warning("SetNpz","Number of points must be >=4 && <= 10000, fNpz set to 4");
661  fNpz = 4;
662  } else if(npz > 10000) {
663  Warning("SetNpz","Number of points must be >=4 && <= 10000, fNpz set to 10000");
664  fNpz = 10000;
665  } else {
666  fNpz = npz;
667  }
668  Update();
669 }
670 
671 ////////////////////////////////////////////////////////////////////////////////
672 /// Initialize the upper and lower bounds to draw the function
673 
675 {
676  fXmin = xmin;
677  fXmax = xmax;
678  fYmin = ymin;
679  fYmax = ymax;
680  fZmin = zmin;
681  fZmax = zmax;
682  Update();
683 }
684 
685 ////////////////////////////////////////////////////////////////////////////////
686 /// Stream an object of class TF3.
687 
688 void TF3::Streamer(TBuffer &R__b)
689 {
690  if (R__b.IsReading()) {
691  UInt_t R__s, R__c;
692  Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
693  if (R__v > 0) {
694  R__b.ReadClassBuffer(TF3::Class(), this, R__v, R__s, R__c);
695  return;
696  }
697 
698  } else {
699  Int_t saved = 0;
700  if (fType != EFType::kFormula && fSave.empty() ) { saved = 1; Save(fXmin,fXmax,fYmin,fYmax,fZmin,fZmax);}
701 
702  R__b.WriteClassBuffer(TF3::Class(),this);
703 
704  if (saved) { fSave.clear(); }
705  }
706 }
707 
708 ////////////////////////////////////////////////////////////////////////////////
709 /// Return x^nx * y^ny * z^nz moment of a 3d function in range [ax,bx],[ay,by],[az,bz]
710 /// \author Gene Van Buren <gene@bnl.gov>
711 
713 {
714  Double_t norm = Integral(ax,bx,ay,by,az,bz,epsilon);
715  if (norm == 0) {
716  Error("Moment3", "Integral zero over range");
717  return 0;
718  }
719 
720  // define integrand function as a lambda : g(x,y,z)= x^(nx) * y^(ny) * z^(nz) * f(x,y,z)
721  auto integrand = [&](double *x, double *) {
722  return std::pow(x[0], nx) * std::pow(x[1], ny) * std::pow(x[2], nz) * this->EvalPar(x, nullptr);
723  };
724  // compute integral of g(x,y,z)
725  TF3 fnc("TF3_ExpValHelper", integrand, ax, bx, ay, by, az, bz, 0);
726  // set same points as current function to get correct max points when computing the integral
727  fnc.fNpx = fNpx;
728  fnc.fNpy = fNpy;
729  fnc.fNpz = fNpz;
730  return fnc.Integral(ax, bx, ay, by, az, bz, epsilon) / norm;
731 }
732 
733 ////////////////////////////////////////////////////////////////////////////////
734 /// Return x^nx * y^ny * z^nz central moment of a 3d function in range [ax,bx],[ay,by],[az,bz]
735 /// \author Gene Van Buren <gene@bnl.gov>
736 
738 {
739  Double_t norm = Integral(ax,bx,ay,by,az,bz,epsilon);
740  if (norm == 0) {
741  Error("CentralMoment3", "Integral zero over range");
742  return 0;
743  }
744 
745  Double_t xbar = 0;
746  Double_t ybar = 0;
747  Double_t zbar = 0;
748  if (nx!=0) {
749  // compute first momentum in x
750  auto integrandX = [&](double *x, double *) { return x[0] * this->EvalPar(x, nullptr); };
751  TF3 fncx("TF3_ExpValHelperx", integrandX, ax, bx, ay, by, az, bz, 0);
752  fncx.fNpx = fNpx;
753  fncx.fNpy = fNpy;
754  fncx.fNpz = fNpz;
755  xbar = fncx.Integral(ax, bx, ay, by, az, bz, epsilon) / norm;
756  }
757  if (ny!=0) {
758  auto integrandY = [&](double *x, double *) { return x[1] * this->EvalPar(x, nullptr); };
759  TF3 fncy("TF3_ExpValHelpery", integrandY, ax, bx, ay, by, az, bz, 0);
760  fncy.fNpx = fNpx;
761  fncy.fNpy = fNpy;
762  fncy.fNpz = fNpz;
763  ybar = fncy.Integral(ax,bx,ay,by,az,bz,epsilon)/norm;
764  }
765  if (nz!=0) {
766  auto integrandZ = [&](double *x, double *) { return x[2] * this->EvalPar(x, nullptr); };
767  TF3 fncz("TF3_ExpValHelperz", integrandZ, ax, bx, ay, by, az, bz, 0);
768  fncz.fNpx = fNpx;
769  fncz.fNpy = fNpy;
770  fncz.fNpz = fNpz;
771  zbar = fncz.Integral(ax,bx,ay,by,az,bz,epsilon)/norm;
772  }
773  // define integrand function as a lambda : g(x,y)= (x-xbar)^(nx) * (y-ybar)^(ny) * f(x,y)
774  auto integrand = [&](double *x, double *) {
775  double xxx = (nx != 0) ? std::pow(x[0] - xbar, nx) : 1.;
776  double yyy = (ny != 0) ? std::pow(x[1] - ybar, ny) : 1.;
777  double zzz = (nz != 0) ? std::pow(x[2] - zbar, nz) : 1.;
778  return xxx * yyy * zzz * this->EvalPar(x, nullptr);
779  };
780  // compute integral of g(x,y, z)
781  TF3 fnc("TF3_ExpValHelper",integrand,ax,bx,ay,by,az,bz,0) ;
782  fnc.fNpx = fNpx;
783  fnc.fNpy = fNpy;
784  fnc.fNpz = fNpz;
785  return fnc.Integral(ax,bx,ay,by,az,bz,epsilon)/norm;
786 }
TF3::SetClippingBoxOn
virtual void SetClippingBoxOn(Double_t xclip=0, Double_t yclip=0, Double_t zclip=0)
Set the function clipping box (for drawing) "on" and define the clipping box.
Definition: TF3.cxx:641
l
auto * l
Definition: textangle.C:4
n
const Int_t n
Definition: legend1.C:16
TF2::Copy
virtual void Copy(TObject &f2) const
Copy this F2 to a new F2.
Definition: TF2.cxx:192
ymax
float ymax
Definition: THbookFile.cxx:95
TF1::ExecuteEvent
virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py)
Execute action corresponding to one event.
Definition: TF1.cxx:1547
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
tt
auto * tt
Definition: textangle.C:16
TVectorD.h
Version_t
short Version_t
Definition: RtypesCore.h:65
TF1::fSave
std::vector< Double_t > fSave
Definition: TF1.h:257
TF1::fHistogram
TH1 * fHistogram
Parent object hooking this function (if one)
Definition: TF1.h:263
f
#define f(i)
Definition: RSha256.hxx:104
Option_t
const char Option_t
Definition: RtypesCore.h:66
TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
TMath::QuietNaN
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition: TMath.h:901
kNPOS
const Ssiz_t kNPOS
Definition: RtypesCore.h:115
TF3::FindMinMax
virtual Double_t FindMinMax(Double_t *x, bool findmax) const
Return minimum/maximum value of the function.
Definition: TF3.cxx:213
TF3::Moment3
virtual Double_t Moment3(Double_t nx, Double_t ax, Double_t bx, Double_t ny, Double_t ay, Double_t by, Double_t nz, Double_t az, Double_t bz, Double_t epsilon=0.000001)
Return x^nx * y^ny * z^nz moment of a 3d function in range [ax,bx],[ay,by],[az,bz].
Definition: TF3.cxx:712
TF3::Copy
virtual void Copy(TObject &f3) const
Copy this F3 to a new F3.
Definition: TF3.cxx:153
TF1::GetParameters
virtual Double_t * GetParameters() const
Definition: TF1.h:519
TString::Data
const char * Data() const
Definition: TString.h:369
TF2::fYmin
Double_t fYmin
Definition: TF2.h:32
TF3::fNpz
Int_t fNpz
Definition: TF3.h:33
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
TNamed::GetTitle
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
TH1::GetPainter
TVirtualHistPainter * GetPainter(Option_t *option="")
Return pointer to painter.
Definition: TH1.cxx:4429
TF1::GetNdim
virtual Int_t GetNdim() const
Definition: TF1.h:484
TBuffer::ReadClassBuffer
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
TF1::GetNpar
virtual Int_t GetNpar() const
Definition: TF1.h:480
r
ROOT::R::TRInterface & r
Definition: Object.C:4
xmax
float xmax
Definition: THbookFile.cxx:95
TObject::Info
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:864
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:890
TF1::fType
EFType fType
Definition: TF1.h:248
TF1::fIntegral
std::vector< Double_t > fIntegral
Definition: TF1.h:258
TVirtualHistPainter.h
TRandom.h
TF3::GetSave
virtual Double_t GetSave(const Double_t *x)
Get value corresponding to X in array of fSave values.
Definition: TF3.cxx:398
Int_t
int Int_t
Definition: RtypesCore.h:45
TF1::DistancetoPrimitive
virtual Int_t DistancetoPrimitive(Int_t px, Int_t py)
Compute distance from point px,py to a function.
Definition: TF1.cxx:1298
TObject::AppendPad
virtual void AppendPad(Option_t *option="")
Append graphics object to current pad.
Definition: TObject.cxx:107
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
x
Double_t x[n]
Definition: legend1.C:17
TF3::TF3
TF3()
coordinates of clipbox
Definition: TF3.cxx:37
TF1::fXmin
Double_t fXmin
Definition: TF1.h:243
TColor.h
TF1::GetMinMaxNDim
virtual Double_t GetMinMaxNDim(Double_t *x, Bool_t findmax, Double_t epsilon=0, Int_t maxiter=0) const
Find the minimum of a function of whatever dimension.
Definition: TF1.cxx:1724
TBuffer
Buffer base class used for serializing objects.
Definition: TBuffer.h:43
TF3::GetRange
virtual void GetRange(Double_t *xmin, Double_t *xmax) const
Return range of a generic N-D function.
Definition: TF1.cxx:2280
TF1::fMethodCall
TMethodCall * fMethodCall
Pointer to histogram used for visualisation.
Definition: TF1.h:264
TColor::SaveColor
static void SaveColor(std::ostream &out, Int_t ci)
Save a color with index > 228 as a C++ statement(s) on output stream out.
Definition: TColor.cxx:2118
TString
Basic string class.
Definition: TString.h:136
TVirtualFitter.h
v
@ v
Definition: rootcling_impl.cxx:3635
b
#define b(i)
Definition: RSha256.hxx:100
bool
TF1::GetParLimits
virtual void GetParLimits(Int_t ipar, Double_t &parmin, Double_t &parmax) const
Return limits for parameter ipar.
Definition: TF1.cxx:1941
TF3::GetMinimumXYZ
virtual Double_t GetMinimumXYZ(Double_t &x, Double_t &y, Double_t &z)
Compute the X, Y and Z values corresponding to the minimum value of the function on its range.
Definition: TF3.cxx:281
TROOT.h
TF3::GetRandom3
virtual void GetRandom3(Double_t &xrandom, Double_t &yrandom, Double_t &zrandom, TRandom *rng=nullptr)
Return 3 random numbers following this function shape.
Definition: TF3.cxx:325
TF3::~TF3
virtual ~TF3()
F3 default destructor.
Definition: TF3.cxx:138
TF1::fNdim
Int_t fNdim
Definition: TF1.h:246
TF3::GetNpz
Int_t GetNpz() const
Definition: TF3.h:94
TF3::SetNpz
virtual void SetNpz(Int_t npz=30)
Set the number of points used to draw the function.
Definition: TF3.cxx:657
TF3::SetClippingBoxOff
virtual void SetClippingBoxOff()
Set the function clipping box (for drawing) "off".
Definition: TF3.cxx:526
TF1::fXmax
Double_t fXmax
Definition: TF1.h:244
TF3::Paint
virtual void Paint(Option_t *option="")
Paint this 3-D function with its current attributes.
Definition: TF3.cxx:501
TBuffer.h
TRandom
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
TF3.h
TF3::fZmin
Double_t fZmin
Definition: TF3.h:31
xmin
float xmin
Definition: THbookFile.cxx:95
TF1::Update
virtual void Update()
Called by functions such as SetRange, SetNpx, SetParameters to force the deletion of the associated h...
Definition: TF1.cxx:3633
h
#define h(i)
Definition: RSha256.hxx:106
TAttLine::GetLineColor
virtual Color_t GetLineColor() const
Return the line color.
Definition: TAttLine.h:33
TF3::Draw
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF3.cxx:177
TF3::fClipBox
Double_t fClipBox[3]
is clip box on
Definition: TF3.h:35
TF1::InitArgs
virtual void InitArgs(const Double_t *x, const Double_t *params)
Initialize parameters addresses.
Definition: TF1.cxx:2471
epsilon
REAL epsilon
Definition: triangle.c:617
a
auto * a
Definition: textangle.C:12
TBuffer::WriteClassBuffer
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
TString::Append
TString & Append(const char *cs)
Definition: TString.h:564
TF3::operator=
TF3 & operator=(const TF3 &rhs)
Operator =.
Definition: TF3.cxx:127
gRandom
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
gDebug
Int_t gDebug
Definition: TROOT.cxx:590
TF3::DistancetoPrimitive
virtual Int_t DistancetoPrimitive(Int_t px, Int_t py)
Compute distance from point px,py to a function.
Definition: TF3.cxx:168
TVirtualPad.h
TF3
A 3-Dim function with parameters.
Definition: TF3.h:28
UInt_t
unsigned int UInt_t
Definition: RtypesCore.h:46
TH1::SetDirectory
virtual void SetDirectory(TDirectory *dir)
By default when an histogram is created, it is added to the list of histogram objects in the current ...
Definition: TH1.cxx:8612
TRandom::Rndm
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
y
Double_t y[n]
Definition: legend1.C:17
TMath::Infinity
Double_t Infinity()
Returns an infinity as defined by the IEEE standard.
Definition: TMath.h:914
TObject::MakeZombie
void MakeZombie()
Definition: TObject.h:49
TF3::CreateHistogram
virtual TH1 * CreateHistogram()
Create a histogram for axis range.
Definition: TF3.cxx:489
TF3::GetMaximumXYZ
virtual Double_t GetMaximumXYZ(Double_t &x, Double_t &y, Double_t &z)
Compute the X, Y and Z values corresponding to the maximum value of the function on its range.
Definition: TF3.cxx:297
TF3::Integral
virtual Double_t Integral(Double_t ax, Double_t bx, Double_t ay, Double_t by, Double_t az, Double_t bz, Double_t epsrel=1.e-6)
Return Integral of a 3d function in range [ax,bx],[ay,by],[az,bz] with a desired relative accuracy.
Definition: TF3.cxx:452
TObject::Warning
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:876
TBuffer::ReadVersion
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
TF1::EvalPar
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=0)
Evaluate function with given coordinates and parameters.
Definition: TF1.cxx:1479
TH3.h
TH3F
3-D histogram with a float per channel (see TH1 documentation)}
Definition: TH3.h:268
ymin
float ymin
Definition: THbookFile.cxx:95
TF2
A 2-Dim function with parameters.
Definition: TF2.h:29
TMath::BinarySearch
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition: TMathBase.h:278
TString::Index
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition: TString.h:639
TVirtualHistPainter::ProcessMessage
virtual void ProcessMessage(const char *mess, const TObject *obj)=0
IntegratorOptions.h
TBuffer::IsReading
Bool_t IsReading() const
Definition: TBuffer.h:86
TF3::fClipBoxOn
Bool_t fClipBoxOn
Definition: TF3.h:34
TF3::fZmax
Double_t fZmax
Definition: TF3.h:32
Double_t
double Double_t
Definition: RtypesCore.h:59
TF1::GetParError
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition: TF1.cxx:1931
TF2::fYmax
Double_t fYmax
Definition: TF2.h:33
TF3::SavePrimitive
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save primitive as a C++ statement(s) on output stream out.
Definition: TF3.cxx:590
TF3::Save
virtual void Save(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Double_t zmin, Double_t zmax)
Save values of function in array fSave.
Definition: TF3.cxx:535
TF3::IsInside
virtual Bool_t IsInside(const Double_t *x) const
Return kTRUE is the point is inside the function range.
Definition: TF3.cxx:478
TObject
Mother of all ROOT objects.
Definition: TObject.h:37
TH1
TH1 is the base class of all histogramm classes in ROOT.
Definition: TH1.h:58
name
char name[80]
Definition: TGX11.cxx:110
ROOT::Math::IntegratorMultiDimOptions::DefaultNCalls
static unsigned int DefaultNCalls()
Definition: IntegratorOptions.cxx:372
gPad
#define gPad
Definition: TVirtualPad.h:287
TAttFill::GetFillColor
virtual Color_t GetFillColor() const
Return the fill area color.
Definition: TAttFill.h:30
TF3::ExecuteEvent
virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py)
Execute action corresponding to one event.
Definition: TF3.cxx:192
TF3::CentralMoment3
virtual Double_t CentralMoment3(Double_t nx, Double_t ax, Double_t bx, Double_t ny, Double_t ay, Double_t by, Double_t nz, Double_t az, Double_t bz, Double_t epsilon=0.000001)
Return x^nx * y^ny * z^nz central moment of a 3d function in range [ax,bx],[ay,by],...
Definition: TF3.cxx:737
TF1::IntegralMultiple
virtual Double_t IntegralMultiple(Int_t n, const Double_t *a, const Double_t *b, Int_t maxpts, Double_t epsrel, Double_t epsabs, Double_t &relerr, Int_t &nfnevl, Int_t &ifail)
This function computes, to an attempted specified accuracy, the value of the integral.
Definition: TF1.cxx:2839
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Class
void Class()
Definition: Class.C:29
TString::ToLower
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
ROOT::Math::ParamFunctorTempl
Param Functor class for Multidimensional functions.
Definition: ParamFunctor.h:273
TF1::GetParameter
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:511
Riostream.h
TF2::fNpy
Int_t fNpy
Definition: TF2.h:34
TF3::Integral
virtual Double_t Integral(Double_t ax, Double_t bx, Double_t ay, Double_t by, Double_t epsrel=1.e-6)
Return Integral of a 2d function in range [ax,bx],[ay,by] with desired relative accuracy (default val...
Definition: TF2.cxx:645
pow
double pow(double, double)
TF1::fNpx
Int_t fNpx
Definition: TF1.h:247
TF1::GetChisquare
Double_t GetChisquare() const
Definition: TF1.h:443
TF3::SetRange
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition: TF3.h:146
TMath::Finite
Int_t Finite(Double_t x)
Check if it is finite with a mask in order to be consistent in presence of fast math.
Definition: TMath.h:771
TH1::Paint
virtual void Paint(Option_t *option="")
Control routine to paint any kind of histograms.
Definition: TH1.cxx:6056
TMath.h
gROOT
#define gROOT
Definition: TROOT.h:406
int