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