Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JRootfit.hh
Go to the documentation of this file.
1 #ifndef __JROOT__JROOTFIT__
2 #define __JROOT__JROOTFIT__
3 
4 #include <vector>
5 #include <set>
6 #include <limits>
7 #include <algorithm>
8 
9 #include "TH1.h"
10 #include "TH2.h"
11 #include "TH3.h"
12 #include "TF1.h"
13 #include "TF2.h"
14 #include "TF3.h"
15 
16 #include "JFit/JGandalf.hh"
17 #include "JMath/JMathSupportkit.hh"
18 #include "JMath/JMathlib.hh"
19 #include "JTools/JRange.hh"
20 
21 
22 /**
23  * \author mdejong
24  */
25 
26 namespace JROOT {}
27 namespace JPP { using namespace JROOT; }
28 
29 namespace JROOT {
30 
31  using JFIT::JGandalf;
33  using JMATH::poisson;
35  using JMATH::getParameter;
37 
38 
39  /**
40  * Type definiton of abscissa range.
41  */
43 
44 
45  /**
46  * Get range of given axis.
47  *
48  * \param pAxis pointer to axis
49  * \return range
50  */
51  inline range_type getRange(TAxis* pAxis)
52  {
53  if (pAxis != NULL)
54  return range_type(pAxis->GetXmin(), pAxis->GetXmax());
55  else
56  return range_type();
57  }
58 
59 
60  /**
61  * Data point for counting measurement.
62  */
63  struct m_count {
64  /**
65  * Default constructor.
66  */
67  m_count() :
68  count(0)
69  {}
70 
71 
72  /**
73  * Constructor.
74  *
75  * \param count count
76  */
77  m_count(const size_t count) :
78  count(count)
79  {}
80 
81 
82  /**
83  * Minimal value for numerical computations.
84  *
85  * \return epsilon
86  */
87  static inline double epsilon()
88  {
89  return std::numeric_limits<double>::min();
90  }
91 
92 
93  /**
94  * Get chi2.
95  *
96  * \param z expectation value
97  * \return chi2
98  */
99  inline double getRho(const double z) const
100  {
101  const double P = poisson(count, z > epsilon() ? z : epsilon());
102 
103  return -log(P > epsilon() ? P : epsilon());
104  }
105 
106 
107  /**
108  * Get derivative of chi2.
109  *
110  * \param z expectation value
111  * \return derivative of chi2
112  */
113  inline double getPsi(const double z) const
114  {
115  return 1.0 - count/(z > epsilon() ? z : epsilon());
116  }
117 
118 
119  size_t count; //!< count
120  };
121 
122 
123  /**
124  * Data point for value with error.
125  */
126  struct m_value {
127  /**
128  * Default constructor.
129  */
131  value(0.0),
132  error(0.0)
133  {}
134 
135 
136  /**
137  * Constructor.
138  *
139  * \param value value
140  * \param error error
141  */
142  m_value(const double value,
143  const double error) :
144  value(value),
145  error (error)
146  {}
147 
148 
149  /**
150  * Get chi2.
151  *
152  * \param z expectation value
153  * \return chi2
154  */
155  inline double getRho(const double z) const
156  {
157  const double u = (value - z) / error;
158 
159  return u*u;
160  }
161 
162 
163  /**
164  * Get derivative of chi2.
165  *
166  * \param z expectation value
167  * \return derivative of chi2
168  */
169  inline double getPsi(const double z) const
170  {
171  const double u = (value - z) / error;
172 
173  return -0.5 * u / error;
174  }
175 
176 
177  double value; //!< value
178  double error; //!< error
179  };
180 
181 
182  /**
183  * 1D data point.
184  */
185  template<class T>
186  struct m_1d :
187  public T
188  {
189  /**
190  * Default constructor.
191  */
192  m_1d() :
193  T(),
194  x(0.0)
195  {}
196 
197 
198  /**
199  * Constructor.
200  *
201  * \param x abscissa
202  * \param v ordinate
203  */
204  m_1d(const double x, const T& v) :
205  T(v),
206  x(x)
207  {}
208 
209  double x; //!< abscissa
210  };
211 
212 
213  /**
214  * 2D data point.
215  */
216  template<class T>
217  struct m_2d :
218  public T
219  {
220  /**
221  * Default constructor.
222  */
223  m_2d() :
224  T(),
225  x(0.0),
226  y(0.0)
227  {}
228 
229 
230  /**
231  * Constructor.
232  *
233  * \param x abscissa
234  * \param y abscissa
235  * \param v ordinate
236  */
237  m_2d(const double x, const double y, const T& v) :
238  T(v),
239  x(x),
240  y(y)
241  {}
242 
243  double x; //!< abscissa
244  double y; //!< abscissa
245  };
246 
247 
248  /**
249  * 3D data point.
250  */
251  template<class T>
252  struct m_3d :
253  public T
254  {
255  /**
256  * Default constructor.
257  */
258  m_3d() :
259  T(),
260  x(0.0),
261  y(0.0),
262  z(0.0)
263  {}
264 
265 
266  /**
267  * Constructor.
268  *
269  * \param x abscissa
270  * \param y abscissa
271  * \param z abscissa
272  * \param v ordinate
273  */
274  m_3d(const double x, const double y, const double z, const T& v) :
275  T(v),
276  x(x),
277  y(y),
278  z(z)
279  {}
280 
281  double x; //!< abscissa
282  double y; //!< abscissa
283  double z; //!< abscissa
284  };
285 
286 
287  /**
288  * Template definition of data structure for set of data points.
289  */
290  template<class T>
291  struct data_type;
292 
293 
294  /**
295  * Template specialisation of data structure for set of 1D data points.
296  */
297  template<class T>
298  struct data_type< m_1d<T> > :
299  public std::vector< m_1d<T> >
300  {
301  /**
302  * Default constructor.
303  */
305  {}
306 
307 
308  /**
309  * Unpack constructor.
310  *
311  * \param h1 1D histogram
312  * \param X abscissa range
313  */
314  data_type(const TH1& h1,
315  const range_type& X = range_type())
316  {
317  unpack(h1, X);
318  }
319 
320 
321  /**
322  * Unpack 1D-histogram.
323  *
324  * \param h1 histogram
325  * \param X abscissa range
326  */
327  void unpack(const TH1& h1,
328  const range_type& X = range_type());
329  };
330 
331 
332  /**
333  * Unpack 1D-histogram.
334  *
335  * \param h1 histogram
336  * \param X abscissa range
337  */
338  template<>
339  void data_type< m_1d<m_count> >::unpack(const TH1& h1,
340  const range_type& X)
341  {
342  for (Int_t ix = 1; ix <= h1.GetXaxis()->GetNbins(); ++ix) {
343 
344  const double x = h1.GetXaxis()->GetBinCenter(ix);
345  const size_t count = h1.GetBinContent(ix);
346 
347  if (X(x)) {
348  this->push_back(m_1d<m_count>(x, count));
349  }
350  }
351  }
352 
353 
354  /**
355  * Unpack 1D-histogram.
356  *
357  * \param h1 histogram
358  * \param X abscissa range
359  */
360  template<>
361  void data_type< m_1d<m_value> >::unpack(const TH1& h1,
362  const range_type& X)
363  {
364  for (Int_t ix = 1; ix <= h1.GetXaxis()->GetNbins(); ++ix) {
365 
366  const double x = h1.GetXaxis()->GetBinCenter(ix);
367  const double value = h1.GetBinContent(ix);
368  const double error = h1.GetBinError (ix);
369 
370  if (X(x) && error > 0.0) {
371  this->push_back(m_1d<m_value>(x, m_value(value, error)));
372  }
373  }
374  }
375 
376 
377  /**
378  * Template specialisation of data structure for set of 2D data points.
379  */
380  template<class T>
381  struct data_type< m_2d<T> > :
382  public std::vector< m_2d<T> >
383  {
384  /**
385  * Default constructor.
386  */
388  {}
389 
390 
391  /**
392  * Unpack constructor.
393  *
394  * \param h2 2D histogram
395  * \param X abscissa range
396  * \param Y abscissa range
397  */
398  data_type(const TH2& h2,
399  const range_type& X = range_type(),
400  const range_type& Y = range_type())
401  {
402  unpack(h2, X, Y);
403  }
404 
405 
406  /**
407  * Unpack 2D-histogram.
408  *
409  * \param h2 histogram
410  * \param X abscissa range
411  * \param Y abscissa range
412  */
413  void unpack(const TH2& h2,
414  const range_type& X = range_type(),
415  const range_type& Y = range_type());
416  };
417 
418 
419  /**
420  * Unpack 2D-histogram.
421  *
422  * \param h2 histogram
423  * \param X abscissa range
424  * \param Y abscissa range
425  */
426  template<>
427  void data_type< m_2d<m_count> >::unpack(const TH2& h2,
428  const range_type& X,
429  const range_type& Y)
430  {
431  for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
432  for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
433 
434  const double x = h2.GetXaxis()->GetBinCenter(ix);
435  const double y = h2.GetYaxis()->GetBinCenter(iy);
436  const size_t count = h2.GetBinContent(ix,iy);
437 
438  if (X(x) && Y(y)) {
439  this->push_back(m_2d<m_count>(x, y, count));
440  }
441  }
442  }
443  }
444 
445 
446  /**
447  * Unpack 2D-histogram.
448  *
449  * \param h2 histogram
450  * \param X abscissa range
451  * \param Y abscissa range
452  */
453  template<>
454  void data_type< m_2d<m_value> >::unpack(const TH2& h2,
455  const range_type& X,
456  const range_type& Y)
457  {
458  for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
459  for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
460 
461  const double x = h2.GetXaxis()->GetBinCenter(ix);
462  const double y = h2.GetYaxis()->GetBinCenter(iy);
463  const double value = h2.GetBinContent(ix,iy);
464  const double error = h2.GetBinError (ix,iy);
465 
466  if (X(x) && Y(y) && error > 0.0) {
467  this->push_back(m_2d<m_value>(x, y, m_value(value, error)));
468  }
469  }
470  }
471  }
472 
473 
474  /**
475  * Template specialisation of data structure for set of 3D data points.
476  */
477  template<class T>
478  struct data_type< m_3d<T> > :
479  public std::vector< m_3d<T> >
480  {
481  /**
482  * Default constructor.
483  */
485  {}
486 
487 
488  /**
489  * Unpack constructor.
490  *
491  * \param h3 2D histogram
492  * \param X abscissa range
493  * \param Y abscissa range
494  * \param Z abscissa range
495  */
496  data_type(const TH3& h3,
497  const range_type& X = range_type(),
498  const range_type& Y = range_type(),
499  const range_type& Z = range_type())
500  {
501  unpack(h3, X, Y, Z);
502  }
503 
504 
505  /**
506  * Unpack 3D-histogram.
507  *
508  * \param h3 histogram
509  * \param X abscissa range
510  * \param Y abscissa range
511  * \param Z abscissa range
512  */
513  void unpack(const TH3& h3,
514  const range_type& X = range_type(),
515  const range_type& Y = range_type(),
516  const range_type& Z = range_type());
517  };
518 
519 
520  /**
521  * Unpack 3D-histogram.
522  *
523  * \param h3 histogram
524  * \param X abscissa range
525  * \param Y abscissa range
526  * \param Z abscissa range
527  */
528  template<>
529  void data_type< m_3d<m_count> >::unpack(const TH3& h3,
530  const range_type& X,
531  const range_type& Y,
532  const range_type& Z)
533  {
534  for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
535  for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
536  for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {
537 
538  const double x = h3.GetXaxis()->GetBinCenter(ix);
539  const double y = h3.GetYaxis()->GetBinCenter(iy);
540  const double z = h3.GetZaxis()->GetBinCenter(iz);
541  const size_t count = h3.GetBinContent(ix,iy,iz);
542 
543  if (X(x) && Y(y) && Z(z)) {
544  this->push_back(m_3d<m_count>(x, y, z, count));
545  }
546  }
547  }
548  }
549  }
550 
551 
552  /**
553  * Unpack 3D-histogram.
554  *
555  * \param h3 histogram
556  * \param X abscissa range
557  * \param Y abscissa range
558  * \param Z abscissa range
559  */
560  template<>
561  void data_type< m_3d<m_value> >::unpack(const TH3& h3,
562  const range_type& X,
563  const range_type& Y,
564  const range_type& Z)
565  {
566  for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
567  for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
568  for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {
569 
570  const double x = h3.GetXaxis()->GetBinCenter(ix);
571  const double y = h3.GetYaxis()->GetBinCenter(iy);
572  const double z = h3.GetZaxis()->GetBinCenter(iz);
573  const double value = h3.GetBinContent(ix,iy,iz);
574  const double error = h3.GetBinError (ix,iy,iz);
575 
576  if (X(x) && Y(y) && Z(z) && error > 0.0) {
577  this->push_back(m_3d<m_value>(x, y, z, m_value(value, error)));
578  }
579  }
580  }
581  }
582  }
583 
584 
585  /**
586  * Wrapper data structure to build ROOT 1D function.
587  */
588  struct JF1 :
589  public TF1
590  {
591  /**
592  * Constructor.
593  *
594  * \param name name
595  * \param f1 function
596  * \param X fit range
597  */
598  template<class JF1_t>
599  JF1(const char* const name,
600  const JF1_t& f1,
601  const range_type& X = range_type()) :
602  TF1(name,
603  helper<JF1_t>(f1),
604  X.getLowerLimit(),
605  X.getUpperLimit(),
606  0)
607  {
608  this->SetNpx(1000);
609  };
610 
611 
612  /**
613  * Helper.
614  */
615  template<class JF1_t>
616  struct helper :
617  public JF1_t
618  {
619  /**
620  * Constructor.
621  *
622  * \param f1 function
623  */
624  helper(const JF1_t& f1) :
625  JF1_t(f1)
626  {}
627 
628 
629  /**
630  * ROOT compatible function call.
631  *
632  * \param x pointer to abscissa values
633  * \param parameters pointer to parameter values
634  * \return function value
635  */
636  double operator()(const double* x, const double* parameters)
637  {
638  setParameters(this, parameters);
639 
640  return this->getValue(x[0]);
641  }
642  };
643  };
644 
645 
646  /**
647  * Wrapper data structure to build ROOT 2D function.
648  */
649  struct JF2 :
650  public TF2
651  {
652  /**
653  * Constructor.
654  *
655  * \param name name
656  * \param f2 function
657  * \param X fit range
658  * \param Y fit range
659  */
660  template<class JF2_t>
661  JF2(const char* const name,
662  const JF2_t& f2,
663  const range_type& X = range_type(),
664  const range_type& Y = range_type()) :
665  TF2(name,
666  helper<JF2_t>(f2),
667  X.getLowerLimit(),
668  X.getUpperLimit(),
669  Y.getLowerLimit(),
670  Y.getUpperLimit(),
671  0)
672  {
673  this->SetNpx(1000);
674  this->SetNpy(1000);
675  };
676 
677 
678  /**
679  * Helper.
680  */
681  template<class JF2_t>
682  struct helper :
683  public JF2_t
684  {
685  /**
686  * Constructor.
687  *
688  * \param f2 function
689  */
690  helper(const JF2_t& f2) :
691  JF2_t(f2)
692  {}
693 
694 
695  /**
696  * ROOT compatible function call.
697  *
698  * \param x pointer to abscissa values
699  * \param parameters pointer to parameter values
700  * \return function value
701  */
702  double operator()(const double* x, const double* parameters)
703  {
704  setParameters(this, parameters);
705 
706  return this->getValue(x[0], x[1]);
707  }
708  };
709  };
710 
711 
712  /**
713  * Wrapper data structure to build ROOT 3D function.
714  */
715  struct JF3 :
716  public TF3
717  {
718  /**
719  * Constructor.
720  *
721  * \param name name
722  * \param f3 function
723  * \param X fit range
724  * \param Y fit range
725  * \param Z fit range
726  */
727  template<class JF3_t>
728  JF3(const char* const name,
729  const JF3_t& f3,
730  const range_type& X = range_type(),
731  const range_type& Y = range_type(),
732  const range_type& Z = range_type()) :
733  TF3(name,
734  helper<JF3_t>(f3),
735  X.getLowerLimit(),
736  X.getUpperLimit(),
737  Y.getLowerLimit(),
738  Y.getUpperLimit(),
739  Z.getLowerLimit(),
740  Z.getUpperLimit(),
741  0)
742  {
743  this->SetNpx(300);
744  this->SetNpy(300);
745  this->SetNpz(300);
746  };
747 
748 
749  /**
750  * Helper.
751  */
752  template<class JF3_t>
753  struct helper :
754  public JF3_t
755  {
756  /**
757  * Constructor.
758  *
759  * \param f3 function
760  */
761  helper(const JF3_t& f3) :
762  JF3_t(f3)
763  {}
764 
765 
766  /**
767  * ROOT compatible function call.
768  *
769  * \param x pointer to abscissa values
770  * \param parameters pointer to parameter values
771  * \return function value
772  */
773  double operator()(const double* x, const double* parameters)
774  {
775  setParameters(this, parameters);
776 
777  return this->getValue(x[0], x[1], x[2]);
778  }
779  };
780  };
781 
782 
783  /**
784  * Auxiliary data structure for list of fixed parameters.
785  */
786  struct index_list :
787  public std::set<size_t>
788  {
789  /**
790  * Default constructor.
791  */
793  {}
794 
795 
796  /**
797  * Constructor.
798  *
799  * \param indices indices
800  */
802  std::set<size_t>(indices)
803  {}
804 
805 
806  /**
807  * Conversion constructor.
808  *
809  * \param parameters parameters
810  */
811  template<class T>
813  {
814  for (size_t i = 0; i != T::parameters.size(); ++i) {
815  if (std::find(parameters.begin(), parameters.end(), T::parameters[i]) != parameters.end()) {
816  this->insert(i);
817  }
818  }
819  }
820  };
821 
822 
823  /**
824  * Base class for result of ROOT Fit.
825  */
826  template<class JFs_t>
827  class JRootfit_t :
828  public JGandalf<JFs_t>
829  {
830  protected:
831  /**
832  * Default constructor.
833  */
835  npx (0),
836  chi2(0.0)
837  {}
838 
839 
840  size_t npx; //!< number of data points
841  double chi2; //!< chi2
842 
843  public:
844  /**
845  * Get function.
846  *
847  * \return function.
848  */
849  const JFs_t& getFunction() const
850  {
851  return this->value;
852  }
853 
854 
855  /**
856  * Get number of parameters.
857  *
858  * \return number of parameters
859  */
860  size_t getNumberOfParameters() const
861  {
862  return JFs_t::parameters.size();
863  }
864 
865 
866  /**
867  * Get number of free parameters.
868  *
869  * \return number of free parameters
870  */
872  {
873  return this->parameters.size();
874  }
875 
876 
877  /**
878  * Get number of data points.
879  *
880  * \return number of data points
881  */
882  size_t getN() const
883  {
884  return npx;
885  }
886 
887 
888  /**
889  * Get chi2.
890  *
891  * \return chi2
892  */
893  double getChi2() const
894  {
895  return chi2;
896  }
897 
898 
899  /**
900  * Get number of degrees of freedom.
901  *
902  * \return number of degrees of freedom
903  */
904  int getNDF() const
905  {
906  return (int) getN() - (int) getNumberOfFreeParameters();
907  }
908 
909 
910  /**
911  * Get value of parameter at given index.
912  *
913  * \param i index
914  */
915  double getValue(size_t i) const
916  {
917  return getParameter(this->value, i);
918  }
919 
920 
921  /**
922  * Get error of parameter at given index.
923  *
924  * \param i index
925  */
926  double getError(size_t i) const
927  {
928  return getParameter(this->error, i);
929  }
930  };
931 
932 
933  /**
934  * ROOT Fit.
935  */
936  template<class JFs_t>
937  class JRootfit :
938  public JRootfit_t<JFs_t>
939  {
940  public:
941 
943 
944 
945  /**
946  * Default constructor.
947  */
949  {}
950 
951 
952  /**
953  * Fit.
954  *
955  * \param h1 histogram
956  * \param f1 start value
957  * \param type type of data for histogram unpacking
958  * \param ls list of fixed parameters
959  * \param X fit range
960  * \return result
961  */
962  template<class T>
963  const result_type& operator()(const TH1& h1,
964  const JFs_t& f1,
965  const T& type,
966  const index_list& ls = index_list(),
967  const range_type& X = range_type())
968  {
969  return eval(f1, ls, data_type< m_1d<T> >(h1, X));
970  }
971 
972 
973  /**
974  * Fit.
975  *
976  * The fitted function is added to the input histogram.
977  *
978  * \param h1 pointer to histogram
979  * \param f1 start value
980  * \param type type of data for histogram unpacking
981  * \param ls list of fixed parameters
982  * \param X fit range
983  * \return result
984  */
985  template<class T>
986  const result_type& operator()(TH1* h1,
987  const JFs_t& f1,
988  const T& type,
989  const index_list& ls = index_list(),
990  const range_type& X = range_type())
991  {
992  (*this)(*h1, f1, type, ls, X);
993 
994  h1->GetListOfFunctions()->Add(new JF1("f1",
995  this->value,
996  JTOOLS::join(X, getRange(h1->GetXaxis()))));
997 
998  return static_cast<const result_type&>(*this);
999  }
1000 
1001 
1002  /**
1003  * Fit.
1004  *
1005  * \param h2 histogram
1006  * \param f2 start value
1007  * \param type type of data for histogram unpacking
1008  * \param ls list of fixed parameters
1009  * \param X fit range
1010  * \param Y fit range
1011  * \return result
1012  */
1013  template<class T>
1014  const result_type& operator()(const TH2& h2,
1015  const JFs_t& f2,
1016  const T& type,
1017  const index_list& ls = index_list(),
1018  const range_type& X = range_type(),
1019  const range_type& Y = range_type())
1020  {
1021  return eval(f2, ls, data_type< m_2d<T> >(h2, X, Y));
1022  }
1023 
1024 
1025  /**
1026  * Fit.
1027  *
1028  * The fitted function is added to the input histogram.
1029  *
1030  * \param h2 pointer to histogram
1031  * \param f2 start value
1032  * \param type type of data for histogram unpacking
1033  * \param ls list of fixed parameters
1034  * \param X fit range
1035  * \param Y fit range
1036  * \return result
1037  */
1038  template<class T>
1039  const result_type& operator()(TH2* h2,
1040  const JFs_t& f2,
1041  const T& type,
1042  const index_list& ls = index_list(),
1043  const range_type& X = range_type(),
1044  const range_type& Y = range_type())
1045  {
1046  (*this)(*h2, f2, type, ls, X, Y);
1047 
1048  h2->GetListOfFunctions()->Add(new JF2("f2",
1049  this->value,
1050  JTOOLS::join(X, getRange(h2->GetXaxis())),
1051  JTOOLS::join(Y, getRange(h2->GetYaxis()))));
1052 
1053  return static_cast<const result_type&>(*this);
1054  }
1055 
1056 
1057  /**
1058  * Fit.
1059  *
1060  * \param h3 histogram
1061  * \param f3 start value
1062  * \param type type of data for histogram unpacking
1063  * \param ls list of fixed parameters
1064  * \param X fit range
1065  * \param Y fit range
1066  * \param Z fit range
1067  * \return result
1068  */
1069  template<class T>
1070  const result_type& operator()(const TH3& h3,
1071  const JFs_t& f3,
1072  const T& type,
1073  const index_list& ls = index_list(),
1074  const range_type& X = range_type(),
1075  const range_type& Y = range_type(),
1076  const range_type& Z = range_type())
1077  {
1078  return eval(f3, ls, data_type< m_3d<T> >(h3, X, Y, Z));
1079  }
1080 
1081 
1082  /**
1083  * Fit.
1084  *
1085  * The fitted function is added to the input histogram.
1086  *
1087  * \param h3 pointer to histogram
1088  * \param f3 start value
1089  * \param type type of data for histogram unpacking
1090  * \param ls list of fixed parameters
1091  * \param X fit range
1092  * \param Y fit range
1093  * \param Z fit range
1094  * \return result
1095  */
1096  template<class T>
1097  const result_type& operator()(TH3* h3,
1098  const JFs_t& f3,
1099  const T& type,
1100  const index_list& ls = index_list(),
1101  const range_type& X = range_type(),
1102  const range_type& Y = range_type(),
1103  const range_type& Z = range_type())
1104  {
1105  (*this)(*h3, f3, type, ls, X, Y, Z);
1106 
1107  h3->GetListOfFunctions()->Add(new JF3("f3",
1108  this->value,
1109  JTOOLS::join(X, getRange(h3->GetXaxis())),
1110  JTOOLS::join(Y, getRange(h3->GetYaxis())),
1111  JTOOLS::join(Z, getRange(h3->GetZaxis()))));
1112 
1113  return static_cast<const result_type&>(*this);
1114  }
1115 
1116 
1117  static JRootfit Fit; //!< Global fit object
1118 
1119  private:
1120  size_t getNumberOfFreeParameters(); // hide method
1121  size_t getN(); // hide method
1122  double getChi2(); // hide method
1123  int getNDF(); // hide method
1124 
1125 
1126  /**
1127  * Evaluate fit.
1128  *
1129  * \param fs start value
1130  * \param ls list of fixed parameters
1131  * \param data data
1132  * \return result
1133  */
1134  template<class T>
1135  const result_type& eval(const JFs_t& fs,
1136  const index_list& ls,
1137  const data_type<T>& data)
1138  {
1139  this->parameters.clear();
1140 
1141  for (size_t i = 0; i != JFs_t::parameters.size(); ++i) {
1142  if (ls.count(i) == 0) {
1143  this->parameters.push_back(JFs_t::parameters[i]);
1144  }
1145  }
1146 
1147  this->value = fs;
1148  this->npx = data.size();
1149  this->chi2 = static_cast<JGandalf<JFs_t>&>(*this)(this->fit, data.begin(), data.end()).chi2;
1150 
1151  return static_cast<const result_type&>(*this);
1152  }
1153 
1154 
1155  /**
1156  * Auxiliary data structure for fit functions.
1157  */
1158  const struct {
1159 
1161 
1162  /**
1163  * Fit function.
1164  *
1165  * \param fs function
1166  * \param mp data point
1167  * \return chi2 and gradient
1168  */
1169  template<class T>
1170  inline const result_type& operator()(const JFs_t& fs, const m_1d<T>& mp) const
1171  {
1172  const double y = fs.getValue(mp.x);
1173 
1174  result.chi2 = mp.getRho(y);
1175  result.gradient = fs.getGradient(mp.x);
1176  result.gradient *= mp.getPsi(y);
1177 
1178  return result;
1179  }
1180 
1181 
1182  /**
1183  * Fit function.
1184  *
1185  * \param fs function
1186  * \param mp data point
1187  * \return chi2 and gradient
1188  */
1189  template<class T>
1190  inline const result_type& operator()(const JFs_t& fs, const m_2d<T>& mp) const
1191  {
1192  const double y = fs.getValue(mp.x, mp.y);
1193 
1194  result.chi2 = mp.getRho(y);
1195  result.gradient = fs.getGradient(mp.x, mp.y);
1196  result.gradient *= mp.getPsi(y);
1197 
1198  return result;
1199  }
1200 
1201 
1202  /**
1203  * Fit function.
1204  *
1205  * \param fs function
1206  * \param mp data point
1207  * \return chi2 and gradient
1208  */
1209  template<class T>
1210  inline const result_type& operator()(const JFs_t& fs, const m_3d<T>& mp) const
1211  {
1212  const double y = fs.getValue(mp.x, mp.y, mp.z);
1213 
1214  result.chi2 = mp.getRho(y);
1215  result.gradient = fs.getGradient(mp.x, mp.y, mp.z);
1216  result.gradient *= mp.getPsi(y);
1217 
1218  return result;
1219  }
1220 
1221  private:
1222  mutable result_type result;
1223  } fit;
1224  };
1225 
1226 
1227  /**
1228  * Global fit object.
1229  */
1230  template<class JFs_t>
1232 
1233 
1234  /**
1235  * Global fit fuction.
1236  *
1237  * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
1238  *
1239  * \param h1 histogram
1240  * \param f1 start value
1241  * \param ls list of fixed parameters
1242  * \param X fit range
1243  * \return result
1244  */
1245  template<class T, class JF1_t>
1246  inline JRootfit_t<JF1_t> Fit(const TH1& h1,
1247  const JF1_t& f1,
1248  const index_list& ls = index_list(),
1249  const range_type& X = range_type())
1250  {
1251  return JRootfit<JF1_t>::Fit(h1, f1, T(), ls, X);
1252  }
1253 
1254 
1255  /**
1256  * Global fit fuction.
1257  *
1258  * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
1259  * The fitted function is added to the input histogram.
1260  *
1261  * \param h1 pointer to histogram
1262  * \param f1 start value
1263  * \param ls list of fixed parameters
1264  * \param X fit range
1265  * \return result
1266  */
1267  template<class T, class JF1_t>
1268  inline JRootfit_t<JF1_t> Fit(TH1* h1,
1269  const JF1_t& f1,
1270  const index_list& ls = index_list(),
1271  const range_type& X = range_type())
1272  {
1273  return JRootfit<JF1_t>::Fit(h1, f1, T(), ls, X);
1274  }
1275 
1276 
1277  /**
1278  * Global fit fuction.
1279  *
1280  * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
1281  *
1282  * \param h2 histogram
1283  * \param f2 start value
1284  * \param ls list of fixed parameters
1285  * \param X fit range
1286  * \param Y fit range
1287  * \return result
1288  */
1289  template<class T, class JF2_t>
1290  inline JRootfit_t<JF2_t> Fit(const TH2& h2,
1291  const JF2_t& f2,
1292  const index_list& ls = index_list(),
1293  const range_type& X = range_type(),
1294  const range_type& Y = range_type())
1295  {
1296  return JRootfit<JF2_t>::Fit(h2, f2, T(), ls, X, Y);
1297  }
1298 
1299 
1300  /**
1301  * Global fit fuction.
1302  *
1303  * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
1304  * The fitted function is added to the input histogram.
1305  *
1306  * \param h2 pointer to histogram
1307  * \param f2 start value
1308  * \param ls list of fixed parameters
1309  * \param X fit range
1310  * \param Y fit range
1311  * \return result
1312  */
1313  template<class T, class JF2_t>
1314  inline JRootfit_t<JF2_t> Fit(TH2* h2,
1315  const JF2_t& f2,
1316  const index_list& ls = index_list(),
1317  const range_type& X = range_type(),
1318  const range_type& Y = range_type())
1319  {
1320  return JRootfit<JF2_t>::Fit(h2, f2, T(), ls, X, Y);
1321  }
1322 
1323 
1324  /**
1325  * Global fit fuction.
1326  *
1327  * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
1328  *
1329  * \param h3 histogram
1330  * \param f3 start value
1331  * \param ls list of fixed parameters
1332  * \param X fit range
1333  * \param Y fit range
1334  * \param Z fit range
1335  * \return result
1336  */
1337  template<class T, class JF3_t>
1338  inline JRootfit_t<JF3_t> Fit(const TH3& h3,
1339  const JF3_t& f3,
1340  const index_list& ls = index_list(),
1341  const range_type& X = range_type(),
1342  const range_type& Y = range_type(),
1343  const range_type& Z = range_type())
1344  {
1345  return JRootfit<JF3_t>::Fit(h3, f3, T(), ls, X, Y, Z);
1346  }
1347 
1348 
1349  /**
1350  * Global fit fuction.
1351  *
1352  * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
1353  * The fitted function is added to the input histogram.
1354  *
1355  * \param h3 pointer to histogram
1356  * \param f3 start value
1357  * \param ls list of fixed parameters
1358  * \param X fit range
1359  * \param Y fit range
1360  * \param Z fit range
1361  * \return result
1362  */
1363  template<class T, class JF3_t>
1364  inline JRootfit_t<JF3_t> Fit(TH3* h3,
1365  const JF3_t& f3,
1366  const index_list& ls = index_list(),
1367  const range_type& X = range_type(),
1368  const range_type& Y = range_type(),
1369  const range_type& Z = range_type())
1370  {
1371  return JRootfit<JF3_t>::Fit(h3, f3, T(), ls, X, Y, Z);
1372  }
1373 }
1374 
1375 #endif
JF3(const char *const name, const JF3_t &f3, const range_type &X=range_type(), const range_type &Y=range_type(), const range_type &Z=range_type())
Constructor.
Definition: JRootfit.hh:728
double getError(size_t i) const
Get error of parameter at given index.
Definition: JRootfit.hh:926
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
int getParameter(const std::string &text)
Get parameter number from text string.
data_type()
Default constructor.
Definition: JRootfit.hh:387
Auxiliary methods for mathematics.
JRange< T, JComparator_t > join(const JRange< T, JComparator_t > &first, const JRange< T, JComparator_t > &second)
Join ranges.
Definition: JRange.hh:659
index_list(const std::initializer_list< double T::* > &parameters)
Conversion constructor.
Definition: JRootfit.hh:812
m_count(const size_t count)
Constructor.
Definition: JRootfit.hh:77
Auxiliary data structure for list of parameters.
Definition: JMathlib.hh:90
JModel_t value
value
Definition: JGandalf.hh:333
const JFs_t & getFunction() const
Get function.
Definition: JRootfit.hh:849
Template definition of data structure for set of data points.
Definition: JRootfit.hh:291
double getRho(const double z) const
Get chi2.
Definition: JRootfit.hh:99
ROOT Fit.
Definition: JRootfit.hh:937
then usage $script< input file >[option] nPossible options count
Definition: JVolume1D.sh:31
const result_type & operator()(const TH2 &h2, const JFs_t &f2, const T &type, const index_list &ls=index_list(), const range_type &X=range_type(), const range_type &Y=range_type())
Fit.
Definition: JRootfit.hh:1014
const result_type & operator()(const TH3 &h3, const JFs_t &f3, const T &type, const index_list &ls=index_list(), const range_type &X=range_type(), const range_type &Y=range_type(), const range_type &Z=range_type())
Fit.
Definition: JRootfit.hh:1070
JRootfit_t()
Default constructor.
Definition: JRootfit.hh:834
m_3d()
Default constructor.
Definition: JRootfit.hh:258
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
m_1d(const double x, const T &v)
Constructor.
Definition: JRootfit.hh:204
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
std::vector< parameter_type > parameters
fit parameters
Definition: JGandalf.hh:330
m_count()
Default constructor.
Definition: JRootfit.hh:67
void unpack(const TH1 &h1, const range_type &X)
Definition: JRootfit.hh:339
const result_type & operator()(const TH1 &h1, const JFs_t &f1, const T &type, const index_list &ls=index_list(), const range_type &X=range_type())
Fit.
Definition: JRootfit.hh:963
data_type(const TH2 &h2, const range_type &X=range_type(), const range_type &Y=range_type())
Unpack constructor.
Definition: JRootfit.hh:398
double operator()(const double *x, const double *parameters)
ROOT compatible function call.
Definition: JRootfit.hh:702
JRootfit()
Default constructor.
Definition: JRootfit.hh:948
double x
abscissa
Definition: JRootfit.hh:243
JRange< int > range_type
double getPsi(const double z) const
Get derivative of chi2.
Definition: JRootfit.hh:169
void setParameters(JF1_t *f1, const double *values)
Set values of all parameters.
Definition: JMathlib.hh:182
double getParameter(const JF1_t &f1, const size_t i)
Get value of parameter at given index.
Definition: JMathlib.hh:169
then set_variable PMT_FILE set_variable DAQ_FILE set_variable OUTPUT_FILE set_variable DETECTOR else fatal Wrong number of arguments fi JPrintTree f $DAQ_FILE type
double z
abscissa
Definition: JRootfit.hh:283
helper(const JF2_t &f2)
Constructor.
Definition: JRootfit.hh:690
data_type(const TH3 &h3, const range_type &X=range_type(), const range_type &Y=range_type(), const range_type &Z=range_type())
Unpack constructor.
Definition: JRootfit.hh:496
helper(const JF1_t &f1)
Constructor.
Definition: JRootfit.hh:624
const JPolynome f1(1.0, 2.0, 3.0)
Function.
2D data point.
Definition: JRootfit.hh:217
Wrapper data structure to build ROOT 1D function.
Definition: JRootfit.hh:588
Auxiliary data structure for list of fixed parameters.
Definition: JRootfit.hh:786
double x
abscissa
Definition: JRootfit.hh:281
Base class for result of ROOT Fit.
Definition: JRootfit.hh:827
static JRootfit Fit
Global fit object.
Definition: JRootfit.hh:1117
data_type(const TH1 &h1, const range_type &X=range_type())
Unpack constructor.
Definition: JRootfit.hh:314
1D data point.
Definition: JRootfit.hh:186
double getRho(const double z) const
Get chi2.
Definition: JRootfit.hh:155
size_t getNumberOfFreeParameters() const
Get number of free parameters.
Definition: JRootfit.hh:871
do set_variable OUTPUT_DIRECTORY $WORKDIR T
JF1(const char *const name, const JF1_t &f1, const range_type &X=range_type())
Constructor.
Definition: JRootfit.hh:599
size_t getN() const
Get number of data points.
Definition: JRootfit.hh:882
size_t getNumberOfParameters() const
Get number of parameters.
Definition: JRootfit.hh:860
JRootfit_t< JF1_t > Fit(const TH1 &h1, const JF1_t &f1, const index_list &ls=index_list(), const range_type &X=range_type())
Global fit fuction.
Definition: JRootfit.hh:1246
const result_type & operator()(TH2 *h2, const JFs_t &f2, const T &type, const index_list &ls=index_list(), const range_type &X=range_type(), const range_type &Y=range_type())
Fit.
Definition: JRootfit.hh:1039
data_type()
Default constructor.
Definition: JRootfit.hh:484
result_type result
Definition: JRootfit.hh:1222
const result_type & operator()(TH3 *h3, const JFs_t &f3, const T &type, const index_list &ls=index_list(), const range_type &X=range_type(), const range_type &Y=range_type(), const range_type &Z=range_type())
Fit.
Definition: JRootfit.hh:1097
m_value()
Default constructor.
Definition: JRootfit.hh:130
then fatal The output file must have the wildcard in the name
Definition: JCanberra.sh:31
range_type getRange(TAxis *pAxis)
Get range of given axis.
Definition: JRootfit.hh:51
m_1d()
Default constructor.
Definition: JRootfit.hh:192
Wrapper data structure to build ROOT 2D function.
Definition: JRootfit.hh:649
Data point for counting measurement.
Definition: JRootfit.hh:63
const result_type & operator()(TH1 *h1, const JFs_t &f1, const T &type, const index_list &ls=index_list(), const range_type &X=range_type())
Fit.
Definition: JRootfit.hh:986
Fit method based on the Levenberg-Marquardt method.
Definition: JGandalf.hh:84
size_t count
count
Definition: JRootfit.hh:119
const result_type & eval(const JFs_t &fs, const index_list &ls, const data_type< T > &data)
Evaluate fit.
Definition: JRootfit.hh:1135
m_2d(const double x, const double y, const T &v)
Constructor.
Definition: JRootfit.hh:237
size_t getNumberOfFreeParameters()
double getChi2()
struct JROOT::JRootfit::@61 fit
Auxiliary data structure for fit functions.
Auxiliary class to define a range between two values.
Wrapper data structure to build ROOT 3D function.
Definition: JRootfit.hh:715
m_value(const double value, const double error)
Constructor.
Definition: JRootfit.hh:142
double x
abscissa
Definition: JRootfit.hh:209
size_t npx
number of data points
Definition: JRootfit.hh:840
double poisson(const size_t n, const double mu)
Poisson probability density distribition.
m_2d()
Default constructor.
Definition: JRootfit.hh:223
static double epsilon()
Minimal value for numerical computations.
Definition: JRootfit.hh:87
data_type()
Default constructor.
Definition: JRootfit.hh:304
double operator()(const double *x, const double *parameters)
ROOT compatible function call.
Definition: JRootfit.hh:773
no fit printf nominal n $STRING awk v X
double y
abscissa
Definition: JRootfit.hh:244
int getNDF() const
Get number of degrees of freedom.
Definition: JRootfit.hh:904
JF2(const char *const name, const JF2_t &f2, const range_type &X=range_type(), const range_type &Y=range_type())
Constructor.
Definition: JRootfit.hh:661
do set_variable MODULE getModule a $WORKDIR detector_a datx L $STRING JEditDetector a $WORKDIR detector_a datx M $MODULE setz o $WORKDIR detector_a datx JEditDetector a $WORKDIR detector_b datx M $MODULE setz o $WORKDIR detector_b datx done echo Output stored at $WORKDIR detector_a datx and $WORKDIR tripod_a txt JDrawDetector2D a $WORKDIR detector_a datx a $WORKDIR detector_b datx L BL o detector $FORMAT $BATCH JDrawDetector2D T $WORKDIR tripod_a txt T $WORKDIR tripod_b txt L BL o tripod $FORMAT $BATCH JCompareDetector a $WORKDIR detector_a datx b $WORKDIR detector_b datx o $WORKDIR abc root &dev null for KEY in X Y Z
double getValue(size_t i) const
Get value of parameter at given index.
Definition: JRootfit.hh:915
Functional algebra.
double f3(const double x, const double y, const double z)
3D function.
Definition: JPolynome3D.cc:23
JModel_t error
error
Definition: JGandalf.hh:334
double operator()(const double *x, const double *parameters)
ROOT compatible function call.
Definition: JRootfit.hh:636
m_3d(const double x, const double y, const double z, const T &v)
Constructor.
Definition: JRootfit.hh:274
double getChi2() const
Get chi2.
Definition: JRootfit.hh:893
helper(const JF3_t &f3)
Constructor.
Definition: JRootfit.hh:761
index_list()
Default constructor.
Definition: JRootfit.hh:792
data_type v[N+1][M+1]
Definition: JPolint.hh:866
double u[N+1]
Definition: JPolint.hh:865
JRootfit_t< JFs_t > result_type
Definition: JRootfit.hh:942
index_list(const std::initializer_list< size_t > &indices)
Constructor.
Definition: JRootfit.hh:801
Data point for value with error.
Definition: JRootfit.hh:126
double y
abscissa
Definition: JRootfit.hh:282
3D data point.
Definition: JRootfit.hh:252
double chi2
chi2
Definition: JRootfit.hh:841
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62
double value
value
Definition: JRootfit.hh:177
JGandalf< JFs_t >::result_type result_type
Definition: JRootfit.hh:1160
double getPsi(const double z) const
Get derivative of chi2.
Definition: JRootfit.hh:113
size_t getNumberOfParameters()
Get number of parameters.
Definition: JMathlib.hh:155
double error
error
Definition: JRootfit.hh:178