Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TFractionFitter.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2// Author: Frank Filthaut F.Filthaut@science.ru.nl 20/05/2002
3// with additions by Bram Wijngaarden <dwijngaa@hef.kun.nl>
4
5/** \class TFractionFitter
6Fits MC fractions to data histogram. À la [HMCMLL](https://cds.cern.ch/record/2296378/files/hbook.pdf),
7see R. Barlow and C. Beeston, Comp. Phys. Comm. 77 (1993) 219-228
8\see https://indico.in2p3.fr/event/2635/contributions/25070/
9\note An alternative interface is described in https://root.cern.ch/doc/master/rf709__BarlowBeeston_8C_source.html
10
11The virtue of this fit is that it takes into account both data and Monte Carlo
12statistical uncertainties. The way in which this is done is through a standard
13likelihood fit using Poisson statistics; however, the template (MC) predictions
14are also varied within statistics, leading to additional contributions to the
15overall likelihood. This leads to many more fit parameters (one per bin per
16template), but the minimisation with respect to these additional parameters is
17done analytically rather than introducing them as formal fit parameters. Some
18special care needs to be taken in the case of bins with zero content. For more
19details please see the original publication cited above.
20
21An example application of this fit is given below. For a TH1* histogram
22("data") fitted as the sum of three Monte Carlo sources ("mc"):
23
24~~~{.cpp}
25{
26 TH1F *data; //data histogram
27 TH1F *mc0; // first MC histogram
28 TH1F *mc1; // second MC histogram
29 TH1F *mc2; // third MC histogram
30 .... // retrieve histograms
31 TObjArray *mc = new TObjArray(3); // MC histograms are put in this array
32 mc->Add(mc0);
33 mc->Add(mc1);
34 mc->Add(mc2);
35 TFractionFitter* fit = new TFractionFitter(data, mc); // initialise
36 fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
37 fit->SetRangeX(1,15); // use only the first 15 bins in the fit
38 Int_t status = fit->Fit(); // perform the fit
39 std::cout << "fit status: " << status << std::endl;
40 if (status == 0) { // check on fit status
41 TH1F* result = (TH1F*) fit->GetPlot();
42 data->Draw("Ep");
43 result->Draw("same");
44 }
45}
46~~~
47
48\note A fully runing example can be found in tutorials/math/fit/fitFraction.C
49
50## Assumptions
51A few assumptions need to be made for the fit procedure to be carried out:
52 1 The total number of events in each template is not too small
53 (so that its Poisson uncertainty can be neglected).
54 2 The number of events in each bin is much smaller than the total
55 number of events in each template (so that multinomial
56 uncertainties can be replaced with Poisson uncertainties).
57
58Biased fit uncertainties may result if these conditions are not fulfilled
59(see e.g. arXiv:0803.2711).
60
61## Instantiation
62A fit object is instantiated through
63 TFractionFitter* fit = new TFractionFitter(data, mc);
64A number of basic checks (intended to ensure that the template histograms
65represent the same "kind" of distribution as the data one) are carried out.
66The TVirtualFitter object is then addressed and all fit parameters (the
67template fractions) declared (initially unbounded).
68
69## Applying constraints
70Fit parameters can be constrained through
71
72 fit->Constrain(parameter #, lower bound, upper bound);
73
74Setting lower bound = upper bound = 0 removes the constraint (a la Minuit);
75however, a function
76
77 fit->Unconstrain(parameter #)
78
79is also provided to simplify this.
80
81## Setting parameter values
82The function
83
84 ROOT::Fit::Fitter* fitter = fit->GetFitter();
85
86is provided for direct access to the ROOT::Fit::Fitter object. This allows to
87set and fix parameter values, limits and set step sizes directly via
88
89 fitter->Config().ParSettings(parameter #).Set(const std::string &name, double value, double step, double lower,
90double upper);
91
92## Restricting the fit range
93The fit range can be restricted through
94
95 fit->SetRangeX(first bin #, last bin #);
96and freed using
97
98 fit->ReleaseRangeX();
99For 2D histograms the Y range can be similarly restricted using
100
101 fit->SetRangeY(first bin #, last bin #);
102 fit->ReleaseRangeY();
103and for 3D histograms also
104
105 fit->SetRangeZ(first bin #, last bin #);
106 fit->ReleaseRangeZ();
107It is also possible to exclude individual bins from the fit through
108
109 fit->ExcludeBin(bin #);
110where the given bin number is assumed to follow the TH1::GetBin() numbering.
111Any bins excluded in this way can be included again using the corresponding
112
113 fit->IncludeBin(bin #);
114
115## Weights histograms
116Weights histograms (for a motivation see the above publication) can be specified
117for the individual MC sources through
118
119 fit->SetWeight(parameter #, pointer to weights histogram);
120and unset by specifying a null pointer.
121
122## Obtaining fit results
123The fit is carried out through
124
125 Int_t status = fit->Fit();
126where status is the code returned from the "MINIMIZE" command. For fits
127that converged, parameter values and errors can be obtained through
128
129 fit->GetResult(parameter #, value, error);
130and the histogram corresponding to the total Monte Carlo prediction (which
131is not the same as a simple weighted sum of the input Monte Carlo distributions)
132can be obtained by
133
134 TH1* result = fit->GetPlot();
135## Using different histograms
136It is possible to change the histogram being fitted through
137
138 fit->SetData(TH1* data);
139and to change the template histogram for a given parameter number through
140
141 fit->SetMC(parameter #, TH1* MC);
142This can speed up code in case of multiple data or template histograms;
143however, it should be done with care as any settings are taken over from
144the previous fit. In addition, neither the dimensionality nor the numbers of
145bins of the histograms should change (in that case it is better to instantiate
146a new TFractionFitter object).
147
148## Errors
149Any serious inconsistency results in an error.
150*/
151
152#include "TH1.h"
153#include "TMath.h"
154#include "TClass.h"
155
156#include "Fit/FitConfig.h"
157#include "Fit/Fitter.h"
158#include "TFitResult.h"
159#include "Math/Functor.h"
160#include "TFractionFitter.h"
161
162
163
164////////////////////////////////////////////////////////////////////////////////
165/// TFractionFitter default constructor.
166
168fFitDone(kFALSE),
169fLowLimitX(0), fHighLimitX(0),
170fLowLimitY(0), fHighLimitY(0),
171fLowLimitZ(0), fHighLimitZ(0),
172fData(nullptr), fIntegralData(0),
173fPlot(nullptr)
174{
175 fFractionFitter = nullptr;
176 fIntegralMCs = nullptr;
177 fFractions = nullptr;
178
179 fNpfits = 0;
180 fNDF = 0;
181 fChisquare = 0;
182 fNpar = 0;
183}
184
185////////////////////////////////////////////////////////////////////////////////
186/// TFractionFitter constructor. Does a complete initialisation (including
187/// consistency checks, default fit range as the whole histogram but without
188/// under- and overflows, and declaration of the fit parameters). Note that
189/// the histograms are not copied, only references are used.
190/// \param[in] data histogram to be fitted
191/// \param[in] MCs array of TH1* corresponding template distributions
192/// \param[in] option can be used to control the print level of the minimization algorithm
193/// - option = "Q" : quite - no message is printed
194/// - option = "V" : verbose - max print out
195/// - option = "" : default: print initial fraction values and result
196
198fFitDone(kFALSE), fChisquare(0), fPlot(nullptr) {
199 fData = data;
200 // Default: include all of the histogram (but without under- and overflows)
201 fLowLimitX = 1;
203 if (fData->GetDimension() > 1) {
204 fLowLimitY = 1;
206 if (fData->GetDimension() > 2) {
207 fLowLimitZ = 1;
209 }
210 }
211 fNpar = MCs->GetEntries();
212 Int_t par;
213 for (par = 0; par < fNpar; ++par) {
214 fMCs.Add(MCs->At(par));
215 // Histogram containing template prediction
216 TString s = TString::Format("Prediction for MC sample %i",par);
217 TH1* pred = (TH1*) ((TH1*)MCs->At(par))->Clone(s);
218 // TFractionFitter manages these histograms
219 pred->SetDirectory(nullptr);
220 pred->SetTitle(s);
221 fAji.Add(pred);
222 }
225
228
230
231 // set print level
232 TString opt(option);
233 opt.ToUpper();
234 if (opt.Contains("Q") ) {
235 fFractionFitter->Config().MinimizerOptions().SetPrintLevel(0);
236 }
237 else if (opt.Contains("V") ) {
238 fFractionFitter->Config().MinimizerOptions().SetPrintLevel(2);
239 }
240 else
241 fFractionFitter->Config().MinimizerOptions().SetPrintLevel(1);
242
244 Double_t defaultStep = 0.01;
245 // set the parameters
246 std::vector<ROOT::Fit::ParameterSettings> & parameters = fFractionFitter->Config().ParamsSettings();
247 parameters.reserve(fNpar);
248 for (par = 0; par < fNpar; ++par) {
249 TString name("frac"); name += par;
250 parameters.push_back(ROOT::Fit::ParameterSettings(name.Data(), defaultFraction, defaultStep) );
251 }
252
253 if (fFractionFitter->Config().MinimizerOptions().ErrorDef() == 1.0 )
254 fFractionFitter->Config().MinimizerOptions().SetErrorDef(0.5);
255
256}
257
258////////////////////////////////////////////////////////////////////////////////
259/// TFractionFitter default destructor
260
263 delete[] fIntegralMCs;
264 delete[] fFractions;
265 if (fPlot) delete fPlot;
266 fAji.Delete();
267}
268
269////////////////////////////////////////////////////////////////////////////////
270/// Change the histogram to be fitted to. Notes:
271/// - Parameter constraints and settings are retained from a possible previous fit.
272/// - Modifying the dimension or number of bins results in an error (in this case
273/// rather instantiate a new TFractionFitter object)
274
280
281////////////////////////////////////////////////////////////////////////////////
282/// Change the histogram for template number `<parm>`. Notes:
283/// - Parameter constraints and settings are retained from a possible previous fit.
284/// - Modifying the dimension or number of bins results in an error (in this case
285/// rather instantiate a new TFractionFitter object)
286
294
295////////////////////////////////////////////////////////////////////////////////
296/// Set bin by bin weights for template number `<parm>` (the parameter numbering
297/// follows that of the input template vector).
298/// Weights can be "unset" by passing a null pointer.
299/// Consistency of the weights histogram with the data histogram is checked at
300/// this point, and an error in case of problems.
301
304 if (fWeights[parm]) {
306 }
307 if (weight) {
308 if (weight->GetNbinsX() != fData->GetNbinsX() ||
309 (fData->GetDimension() > 1 && weight->GetNbinsY() != fData->GetNbinsY()) ||
310 (fData->GetDimension() > 2 && weight->GetNbinsZ() != fData->GetNbinsZ())) {
311 Error("SetWeight","Inconsistent weights histogram for source %d", parm);
312 return;
313 }
314 TString ts = "weight hist: "; ts += weight->GetName();
315 fWeights.AddAt(weight,parm);
316 }
317}
318
319////////////////////////////////////////////////////////////////////////////////
320/// Give direct access to the underlying fitter class. This can be
321/// used e.g. to modify parameter values or step sizes.
322
326
327////////////////////////////////////////////////////////////////////////////////
328/// Function for internal use, checking parameter validity
329/// An invalid parameter results in an error.
330
333 Error("CheckParNo","Invalid parameter number %d",parm);
334 }
335}
336
337////////////////////////////////////////////////////////////////////////////////
338/// Set the X range of the histogram to be used in the fit.
339/// Use ReleaseRangeX() to go back to fitting the full histogram.
340/// The consistency check ensures that no empty fit range occurs (and also
341/// recomputes the bin content integrals).
342/// \param[in] low lower X bin number
343/// \param[in] high upper X bin number
344
346 fLowLimitX = (low > 0 ) ? low : 1;
347 fHighLimitX = ( high > 0 && high <= fData->GetNbinsX()) ? high : fData->GetNbinsX();
349}
350
351////////////////////////////////////////////////////////////////////////////////
352/// Release restrictions on the X range of the histogram to be used in the fit.
353
359
360////////////////////////////////////////////////////////////////////////////////
361/// Set the Y range of the histogram to be used in the fit (2D or 3D histograms only).
362/// Use ReleaseRangeY() to go back to fitting the full histogram.
363/// The consistency check ensures that no empty fit range occurs (and also
364/// recomputes the bin content integrals).
365/// \param[in] low lower X bin number
366/// \param[in] high upper X bin number
367
369 if (fData->GetDimension() < 2) {
370 Error("SetRangeY","Y range cannot be set for 1D histogram");
371 return;
372 }
373
374 fLowLimitY = (low > 0) ? low : 1;
375 fHighLimitY = (high > 0 && high <= fData->GetNbinsY()) ? high : fData->GetNbinsY();
377}
378
379////////////////////////////////////////////////////////////////////////////////
380/// Release restrictions on the Y range of the histogram to be used in the fit.
381
387
388
389////////////////////////////////////////////////////////////////////////////////
390/// Set the Z range of the histogram to be used in the fit (3D histograms only).
391/// Use ReleaseRangeY() to go back to fitting the full histogram.
392/// The consistency check ensures that no empty fit range occurs (and also
393/// recomputes the bin content integrals).
394/// \param[in] low lower X bin number
395/// \param[in] high upper X bin number
396
398 if (fData->GetDimension() < 3) {
399 Error("SetRangeZ","Z range cannot be set for 1D or 2D histogram");
400 return;
401 }
402
403
404 fLowLimitZ = (low > 0) ? low : 1;
405 fHighLimitZ = (high > 0 && high <= fData->GetNbinsZ()) ? high : fData->GetNbinsZ();
407}
408
409////////////////////////////////////////////////////////////////////////////////
410/// Release restrictions on the Z range of the histogram to be used in the fit.
411
417
418////////////////////////////////////////////////////////////////////////////////
419/// Exclude the given bin from the fit. The bin numbering to be used is that
420/// of TH1::GetBin().
421
423 int excluded = fExcludedBins.size();
424 for (int b = 0; b < excluded; ++b) {
425 if (fExcludedBins[b] == bin) {
426 Error("ExcludeBin", "bin %d already excluded", bin);
427 return;
428 }
429 }
430 fExcludedBins.push_back(bin);
431 // This call serves to properly (re)determine the number of degrees of freeom
433}
434
435////////////////////////////////////////////////////////////////////////////////
436/// Include the given bin in the fit, if it was excluded before using ExcludeBin().
437/// The bin numbering to be used is that of TH1::GetBin().
438
440 for (std::vector<Int_t>::iterator it = fExcludedBins.begin();
441 it != fExcludedBins.end(); ++it) {
442 if (*it == bin) {
443 fExcludedBins.erase(it);
444 // This call serves to properly (re)determine the number of degrees of freeom
446 return;
447 }
448 }
449 Error("IncludeBin", "bin %d was not excluded", bin);
450}
451
452////////////////////////////////////////////////////////////////////////////////
453/// Function for internal use, checking whether the given bin is
454/// excluded from the fit or not.
455
457 for (unsigned int b = 0; b < fExcludedBins.size(); ++b)
458 if (fExcludedBins[b] == bin) return true;
459 return false;
460}
461
462////////////////////////////////////////////////////////////////////////////////
463/// Constrain the values of parameter number `<parm>` (the parameter numbering
464/// follows that of the input template vector).
465/// Use UnConstrain() to remove this constraint.
466
469 assert( parm >= 0 && parm < (int) fFractionFitter->Config().ParamsSettings().size() );
470 fFractionFitter->Config().ParSettings(parm).SetLimits(low,high);
471}
472
473////////////////////////////////////////////////////////////////////////////////
474/// Remove the constraints on the possible values of parameter `<parm>`.
475
478 fFractionFitter->Config().ParSettings(parm).RemoveLimits();
479}
480
481////////////////////////////////////////////////////////////////////////////////
482/// Function used internally to check the consistency between the
483/// various histograms. Checks are performed on nonexistent or empty
484/// histograms, the precise histogram class, and the number of bins.
485/// In addition, integrals over the "allowed" bin ranges are computed.
486/// Any inconsistency results in a error.
487
489 if (! fData) {
490 Error("CheckConsistency","Nonexistent data histogram");
491 return;
492 }
494 Int_t x,y,z,par;
496 fIntegralData = 0;
497 fNpfits = 0;
498 for (z = minZ; z <= maxZ; ++z) {
499 for (y = minY; y <= maxY; ++y) {
500 for (x = minX; x <= maxX; ++x) {
501 if (IsExcluded(fData->GetBin(x, y, z))) continue;
502 fNpfits++;
504 }
505 }
506 }
507 if (fIntegralData <= 0) {
508 Error("CheckConsistency","Empty data histogram");
509 return;
510 }
511 TClass* cl = fData->Class();
512
513 fNDF = fNpfits - fNpar;
514
515 if (fNpar < 2) {
516 Error("CheckConsistency","Need at least two MC histograms");
517 return;
518 }
519
520 for (par = 0; par < fNpar; ++par) {
521 TH1 *h = (TH1*)fMCs.At(par);
522 if (! h) {
523 Error("CheckConsistency","Nonexistent MC histogram for source #%d",par);
524 return;
525 }
526 if ((! h->Class()->InheritsFrom(cl)) || h->GetNbinsX() != fData->GetNbinsX() ||
527 (fData->GetDimension() > 1 && h->GetNbinsY() != fData->GetNbinsY()) ||
528 (fData->GetDimension() > 2 && h->GetNbinsZ() != fData->GetNbinsZ())) {
529 Error("CheckConsistency","Histogram inconsistency for source #%d",par);
530 return;
531 }
532 fIntegralMCs[par] = 0;
533 for (z = minZ; z <= maxZ; ++z) {
534 for (y = minY; y <= maxY; ++y) {
535 for (x = minX; x <= maxX; ++x) {
536 Int_t bin = fData->GetBin(x, y, z);
537 if (IsExcluded(bin)) continue;
538 Double_t MCEvents = h->GetBinContent(bin);
539 if (MCEvents < 0) {
540 Error("CheckConsistency", "Number of MC events (bin = %d, par = %d) cannot be negative: "
541 " their distribution is binomial (see paper)", bin, par);
542 }
543 fIntegralMCs[par] += MCEvents;
544 }
545 }
546 }
547 if (fIntegralMCs[par] <= 0) {
548 Error("CheckConsistency","Empty MC histogram #%d",par);
549 }
550 }
551}
552
553////////////////////////////////////////////////////////////////////////////////
554/// Perform the fit with the default UP value.
555/// The value returned is the minimisation status.
556
558
559 // remove any existing output histogram
560 if (fPlot) {
561 delete fPlot; fPlot = nullptr;
562 }
563
564 // Make sure the correct likelihood computation is used
567
568 // fit
569 Bool_t status = fFractionFitter->FitFCN();
570 if (!status) Warning("Fit","Abnormal termination of minimization.");
571
572 fFitDone = kTRUE;
573
574 // determine goodness of fit
576
577 // create a new result class
579 TString name = TString::Format("TFractionFitter_result_of_%s",fData->GetName() );
580 fr->SetName(name); fr->SetTitle(name);
581 return TFitResultPtr(fr);
582}
583
584////////////////////////////////////////////////////////////////////////////////
585/// Set UP to the given value (see class TMinuit), and perform a MINOS minimisation.
586
588 if (! fFitDone) {
589 Error("ErrorAnalysis","Fit not yet performed");
590 return;
591 }
592
593
594 Double_t up = UP > 0 ? UP : 0.5;
595 fFractionFitter->Config().MinimizerOptions().SetErrorDef(up);
597 if (!status) {
598 Error("ErrorAnalysis","Error return from MINOS: %d",fFractionFitter->Result().Status());
599 }
600}
601
602////////////////////////////////////////////////////////////////////////////////
603/// Obtain the fit result for parameter `<parm>` (the parameter numbering
604/// follows that of the input template vector).
605
608 if (! fFitDone) {
609 Error("GetResult","Fit not yet performed");
610 return;
611 }
612 value = fFractionFitter->Result().Parameter(parm);
613 error = fFractionFitter->Result().Error(parm);
614}
615
616////////////////////////////////////////////////////////////////////////////////
617/// Return the "template prediction" corresponding to the fit result (this is not
618/// the same as the weighted sum of template distributions, as template statistical
619/// uncertainties are taken into account).
620/// Note that the name of this histogram will simply be the same as that of the
621/// "data" histogram, prefixed with the string "Fraction fit to hist: ".
622/// Note also that the histogram is managed by the TFractionFitter class, so the returned pointer will be invalid if
623/// the class is deleted
624
626 if (! fFitDone) {
627 Error("GetPlot","Fit not yet performed");
628 return nullptr;
629 }
630 if (! fPlot) {
631 Double_t f = 0;
632 const Double_t * par = fFractionFitter->Result().GetParams();
633 assert(par);
634 ComputeFCN(f, par, 3);
635 }
636 return fPlot;
637}
638
639////////////////////////////////////////////////////////////////////////////////
640/// Used internally to obtain the bin ranges according to the dimensionality of
641/// the histogram and the limits set by hand.
642
644 Int_t& minZ, Int_t& maxZ) const {
645 if (fData->GetDimension() < 2) {
646 minY = maxY = minZ = maxZ = 0;
649 } else if (fData->GetDimension() < 3) {
650 minZ = maxZ = 0;
655 } else {
662 }
663}
664
665////////////////////////////////////////////////////////////////////////////////
666/// Used internally to compute the likelihood value.
667
669{
670 // normalise the fit parameters
671 Int_t bin, mc;
673 Int_t x,y,z;
675 for (mc = 0; mc < fNpar; ++mc) {
677 TH1 *h = (TH1*)fMCs[mc];
678 TH1 *hw = (TH1*)fWeights[mc];
679 if (hw) {
680 tot = 0;
681 for (z = minZ; z <= maxZ; ++z) {
682 for (y = minY; y <= maxY; ++y) {
683 for (x = minX; x <= maxX; ++x) {
684 if (IsExcluded(fData->GetBin(x, y, z))) continue;
685 Double_t weight = hw->GetBinContent(x, y, z);
686 if (weight <= 0) {
687 Error("ComputeFCN","Invalid weight encountered for MC source %d",mc);
688 return;
689 }
690 tot += weight * h->GetBinContent(x, y, z);
691 }
692 }
693 }
694 } else tot = fIntegralMCs[mc];
696 }
697
698 if (flag == 3) {
699 TString ts = "Fraction fit to hist: "; ts += fData->GetName();
700 fPlot = (TH1*) fData->Clone(ts.Data());
701 // plot histogram is managed by TFractionFitter
702 fPlot->SetDirectory(nullptr);
703 fPlot->Reset();
704 }
705 // likelihood computation
706 Double_t result = 0;
707 for (z = minZ; z <= maxZ; ++z) {
708 for (y = minY; y <= maxY; ++y) {
709 for (x = minX; x <= maxX; ++x) {
710 bin = fData->GetBin(x, y, z);
711 if (IsExcluded(bin)) continue;
712
713 // Solve for the "predictions"
714 int k0 = 0;
715 Double_t ti = 0.0; Double_t aki = 0.0;
716 FindPrediction(bin, ti, k0, aki);
717
719 for (mc = 0; mc < fNpar; ++mc) {
720 TH1 *h = (TH1*)fMCs[mc];
721 TH1 *hw = (TH1*)fWeights[mc];
723 Double_t binContent = h->GetBinContent(bin);
724 Double_t weight = hw ? hw->GetBinContent(bin) : 1;
725 if (k0 >= 0 && fFractions[mc] == fFractions[k0]) {
727 } else {
728 binPrediction = binContent > 0 ? binContent / (1+weight*fFractions[mc]*ti) : 0;
729 }
730
733 if (binContent > 0 && binPrediction > 0)
735
736 if (flag == 3) {
737 ((TH1*)fAji.At(mc))->SetBinContent(bin, binPrediction);
738 }
739 }
740
741 if (flag == 3) {
743 }
744
746 Double_t found = fData->GetBinContent(bin);
747 if (found > 0 && prediction > 0)
748 result += found*TMath::Log(prediction);
749 }
750 }
751 }
752
753 f = -result;
754}
755
756////////////////////////////////////////////////////////////////////////////////
757/// Function used internally to obtain the template prediction in the individual bins
758/// 'bin' <=> 'i' (paper)
759/// 'par' <=> 'j' (paper)
760
762 std::vector<Double_t> wgtFrac(fNpar); // weighted fractions (strengths of the sources)
763 std::vector<Double_t> a_ji(fNpar); // number of observed MC events for bin 'i' and par (source) 'j'
764 Double_t d_i = fData->GetBinContent(bin); // number of events in the real data for bin 'i'
765
766 // Cache the weighted fractions and the number of observed MC events
767 // Sanity check: none of the fractions should be == 0
768 for (Int_t par = 0; par < fNpar; ++par) {
769 a_ji[par] = ((TH1*)fMCs.At(par))->GetBinContent(bin);
770 TH1* hw = (TH1*)fWeights.At(par);
771 wgtFrac[par] = hw ? hw->GetBinContent(bin) * fFractions[par] : fFractions[par];
772 if (wgtFrac[par] == 0) {
773 Error("FindPrediction", "Fraction[%d] = 0!", par);
774 return;
775 }
776 }
777
778 // Case data = 0
779 if (TMath::Nint(d_i) == 0) {
780 t_i = 1;
781 k_0 = -1;
782 A_ki = 0;
783 return;
784 }
785
786 // Case one or more of the MC bin contents == 0 -> find largest fraction
787 // k_0 stores the source index of the largest fraction
788 k_0 = 0;
790 for (Int_t par = 1; par < fNpar; ++par) {
791 if (wgtFrac[par] > maxWgtFrac) {
792 k_0 = par;
793 maxWgtFrac = wgtFrac[par];
794 }
795 }
796 Double_t t_min = -1 / maxWgtFrac; // t_i cannot be smaller than this value (see paper, par 5)
797
798 // Determine if there are more sources which have the same maximum contribution (fraction)
800 for (Int_t par = 0; par < fNpar; ++par) {
801 if (par == k_0) continue;
802 if (wgtFrac[par] == maxWgtFrac) {
803 nMax++;
804 contentsMax += a_ji[par];
805 }
806 }
807
808 // special action if there is a zero in the number of entries for the MC source with
809 // the largest strength (fraction) -> See Paper, Paragraph 5
810 if (contentsMax == 0) {
811 A_ki = d_i / (1.0 + maxWgtFrac);
812 for (Int_t par = 0; par < fNpar; ++par) {
813 if (par == k_0 || wgtFrac[par] == maxWgtFrac) continue;
814 A_ki -= a_ji[par] * wgtFrac[par] / (maxWgtFrac - wgtFrac[par]);
815 }
816 if (A_ki > 0) {
817 A_ki /= nMax;
818 t_i = t_min;
819 return;
820 }
821 }
822 k_0 = -1;
823
824 // Case of nonzero histogram contents: solve for t_i using Newton's method
825 // The equation that needs to be solved:
826 // func(t_i) = \sum\limits_j{\frac{ p_j a_{ji} }{1 + p_j t_i}} - \frac{d_i}{1 - t_i} = 0
827 t_i = 0; Double_t step = 0.2;
828 Int_t maxIter = 100000; // maximum number of iterations
829 for(Int_t i = 0; i < maxIter; ++i) {
830 if (t_i >= 1 || t_i < t_min) {
831 step /= 10;
832 t_i = 0;
833 }
834 Double_t func = - d_i / (1.0 - t_i);
835 Double_t deriv = func / (1.0 - t_i);
836 for (Int_t par = 0; par < fNpar; ++par) {
837 Double_t r = 1.0 / (t_i + 1.0 / wgtFrac[par]);
838 func += a_ji[par] * r;
839 deriv -= a_ji[par] * r * r;
840 }
841 if (TMath::Abs(func) < 1e-12) return; // solution found
842 Double_t delta = - func / deriv; // update delta
843 if (TMath::Abs(delta) > step)
844 delta = (delta > 0) ? step : -step; // correct delta if it becomes too large
845 t_i += delta;
846 if (TMath::Abs(delta) < 1e-13) return; // solution found
847 } // the loop breaks when the solution is found, or when the maximum number of iterations is exhausted
848
849 Warning("FindPrediction", "Did not find solution for t_i in %d iterations", maxIter);
850
851 return;
852}
853
854#ifdef OLD
855////////////////////////////////////////////////////////////////////////////////
856/// Function called by the minimisation package. The actual functionality is passed
857/// on to the TFractionFitter::ComputeFCN member function.
858
860 TFractionFitter* fitter = dynamic_cast<TFractionFitter*>(fFractionFitter->GetObjectFit());
861 if (!fitter) {
862 Error("TFractionFitFCN","Invalid fit object encountered!");
863 return;
864 }
865 fitter->ComputeFCN(npar, gin, f, par, flag);
866}
867#endif
868
869////////////////////////////////////////////////////////////////////////////////
870/// Return the likelihood ratio Chi-squared (chi2) for the fit.
871/// The value is computed when the fit is executed successfully.
872/// Chi2 calculation is based on the "likelihood ratio" lambda,
873/// lambda = L(y;n) / L(m;n),
874/// where L(y;n) is the likelihood of the fit result `<y>` describing the data `<n>`
875/// and L(m;n) is the likelihood of an unknown "true" underlying distribution
876/// `<m>` describing the data `<n>`. Since `<m>` is unknown, the data distribution is
877/// used instead,
878/// lambda = L(y;n) / L(n;n).
879/// Note that this ratio is 1 if the fit is perfect. The chi2 value is then
880/// computed according to
881/// chi2 = -2*ln(lambda).
882/// This parameter can be shown to follow a Chi-square distribution. See for
883/// example S. Baker and R. Cousins, "Clarification of the use of chi-square
884/// and likelihood functions in fits to histograms", Nucl. Instr. Meth. A221,
885/// pp. 437-442 (1984)
886
891
892////////////////////////////////////////////////////////////////////////////////
893/// return the number of degrees of freedom in the fit
894/// the fNDF parameter has been previously computed during a fit.
895/// The number of degrees of freedom corresponds to the number of points
896/// used in the fit minus the number of templates.
897
899{
900 if (fNDF == 0) return fNpfits-fNpar;
901 return fNDF;
902}
903
904////////////////////////////////////////////////////////////////////////////////
905/// return the fit probability
906
908{
910 if (ndf <= 0) return 0;
911 return TMath::Prob(fChisquare,ndf);
912}
913
914////////////////////////////////////////////////////////////////////////////////
915/// Method used internally to compute the likelihood ratio chi2
916/// See the function GetChisquare() for details
917
919{
920 if ( !fFitDone ) {
921 Error("ComputeChisquareLambda","Fit not yet (successfully) performed");
922 fChisquare = 0;
923 return;
924 }
925
926 // fPlot must be initialized and filled. Leave this to the GetPlot() method.
927 if (! fPlot)
928 GetPlot();
929
932
933 Double_t logLyn = 0; // likelihood of prediction
934 Double_t logLmn = 0; // likelihood of data ("true" distribution)
935 for(Int_t x = minX; x <= maxX; x++) {
936 for(Int_t y = minY; y <= maxY; y++) {
937 for(Int_t z = minZ; z <= maxZ; z++) {
938 if (IsExcluded(fData->GetBin(x, y, z))) continue;
941 if(fi != 0) logLyn += di * TMath::Log(fi) - fi;
942 if(di != 0) logLmn += di * TMath::Log(di) - di;
943 for(Int_t j = 0; j < fNpar; j++) {
944 Double_t aji = ((TH1*)fMCs.At(j))->GetBinContent(x, y, z);
945 Double_t bji = ((TH1*)fAji.At(j))->GetBinContent(x, y, z);
946 if(bji != 0) logLyn += aji * TMath::Log(bji) - bji;
947 if(aji != 0) logLmn += aji * TMath::Log(aji) - aji;
948 }
949 }
950 }
951 }
952
953 fChisquare = -2*logLyn + 2*logLmn;
954
955 return;
956}
957
958////////////////////////////////////////////////////////////////////////////////
959/// Return the adjusted MC template (Aji) for template (parm).
960/// Note that the (Aji) times fractions only sum to the total prediction
961/// of the fit if all weights are 1.
962/// Note also that the histogram is managed by the TFractionFitter class, so the returned pointer will be invalid if
963/// the class is deleted
964
966{
968 if ( !fFitDone ) {
969 Error("GetMCPrediction","Fit not yet performed");
970 return nullptr;
971 }
972 return (TH1*) fAji.At(parm);
973}
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
void TFractionFitFCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag)
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
char name[80]
Definition TGX11.cxx:110
Double_t GetBinContent(Int_t bin) const override
Fitter class, entry point for performing all type of fits.
Definition Fitter.h:79
const FitResult & Result() const
get fit result
Definition Fitter.h:397
bool FitFCN(unsigned int npar, Function &fcn, const double *params=nullptr, unsigned int dataSize=0, int fitType=0)
Fit using the a generic FCN function as a C++ callable object implementing double () (const double *)...
Definition Fitter.h:647
const FitConfig & Config() const
access to the fit configuration (const method)
Definition Fitter.h:425
bool SetFCN(unsigned int npar, Function &fcn, const double *params=nullptr, unsigned int dataSize=0, int fitType=0)
Set a generic FCN function as a C++ callable object implementing double () (const double *) Note that...
Definition Fitter.h:654
bool CalculateMinosErrors()
perform an error analysis on the result using MINOS To be called only after fitting and when a minimi...
Definition Fitter.cxx:590
Class, describing value, limits and step size of the parameters Provides functionality also to set/re...
Documentation for class Functor class.
Definition Functor.h:49
Documentation for the abstract class IBaseFunctionMultiDim.
Definition IFunction.h:63
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition TClass.h:84
Provides an indirection to the TFitResult class and with a semantics identical to a TFitResult pointe...
Extends the ROOT::Fit::Result class with a TNamed inheritance providing easy possibility for I/O.
Definition TFitResult.h:34
Fits MC fractions to data histogram.
void ComputeFCN(Double_t &f, const Double_t *par, Int_t flag)
Used internally to compute the likelihood value.
void GetRanges(Int_t &minX, Int_t &maxX, Int_t &minY, Int_t &maxY, Int_t &minZ, Int_t &maxZ) const
Used internally to obtain the bin ranges according to the dimensionality of the histogram and the lim...
void FindPrediction(int bin, double &t_i, int &k_0, double &A_ki) const
Function used internally to obtain the template prediction in the individual bins 'bin' <=> 'i' (pape...
TH1 * GetPlot()
Return the "template prediction" corresponding to the fit result (this is not the same as the weighte...
Int_t fLowLimitX
First bin in X dimension.
Double_t GetProb() const
return the fit probability
void GetResult(Int_t parm, Double_t &value, Double_t &error) const
Obtain the fit result for parameter <parm> (the parameter numbering follows that of the input templat...
TObjArray fAji
Array of pointers to predictions of real template distributions.
void CheckConsistency()
Function used internally to check the consistency between the various histograms.
Bool_t fFitDone
Flags whether a valid fit has been performed.
void SetRangeX(Int_t low, Int_t high)
Set the X range of the histogram to be used in the fit.
void SetMC(Int_t parm, TH1 *MC)
Change the histogram for template number <parm>.
TObjArray fMCs
Array of pointers to template histograms.
Double_t EvaluateFCN(const Double_t *par)
ROOT::Fit::Fitter * GetFitter() const
Give direct access to the underlying fitter class.
void SetRangeY(Int_t low, Int_t high)
Set the Y range of the histogram to be used in the fit (2D or 3D histograms only).
~TFractionFitter() override
TFractionFitter default destructor.
void ExcludeBin(Int_t bin)
Exclude the given bin from the fit.
TFractionFitter()
TFractionFitter default constructor.
TObjArray fWeights
Array of pointers to corresponding weight factors (may be null)
void IncludeBin(Int_t bin)
Include the given bin in the fit, if it was excluded before using ExcludeBin().
ROOT::Fit::Fitter * fFractionFitter
Pointer to Fitter class.
void Constrain(Int_t parm, Double_t low, Double_t high)
Constrain the values of parameter number <parm> (the parameter numbering follows that of the input te...
void SetData(TH1 *data)
Change the histogram to be fitted to.
TH1 * fPlot
Pointer to histogram containing summed template predictions.
Int_t fHighLimitY
Last bin in Y dimension.
bool IsExcluded(Int_t bin) const
Function for internal use, checking whether the given bin is excluded from the fit or not.
void SetRangeZ(Int_t low, Int_t high)
Set the Z range of the histogram to be used in the fit (3D histograms only).
Int_t fNpar
number of fit parameters
Int_t fHighLimitZ
Last bin in Z dimension.
Int_t fNpfits
Number of points used in the fit.
Double_t GetChisquare() const
Return the likelihood ratio Chi-squared (chi2) for the fit.
void UnConstrain(Int_t parm)
Remove the constraints on the possible values of parameter <parm>.
Int_t fLowLimitY
First bin in Y dimension.
Int_t fHighLimitX
Last bin in X dimension.
Int_t GetNDF() const
return the number of degrees of freedom in the fit the fNDF parameter has been previously computed du...
Double_t * fFractions
Template fractions scaled to the "data" histogram statistics.
Int_t fLowLimitZ
First bin in Z dimension.
void ReleaseRangeZ()
Release restrictions on the Z range of the histogram to be used in the fit.
void SetWeight(Int_t parm, TH1 *weight)
Set bin by bin weights for template number <parm> (the parameter numbering follows that of the input ...
Double_t fIntegralData
"data" histogram content integral over allowed fit range
void ReleaseRangeX()
Release restrictions on the X range of the histogram to be used in the fit.
void ErrorAnalysis(Double_t UP)
Set UP to the given value (see class TMinuit), and perform a MINOS minimisation.
void ComputeChisquareLambda()
Method used internally to compute the likelihood ratio chi2 See the function GetChisquare() for detai...
Int_t fNDF
Number of degrees of freedom in the fit.
TFitResultPtr Fit()
Perform the fit with the default UP value.
Double_t * fIntegralMCs
Same for template histograms (weights not taken into account)
Double_t fChisquare
Template fit chisquare.
void ReleaseRangeY()
Release restrictions on the Y range of the histogram to be used in the fit.
TH1 * fData
Pointer to the "data" histogram to be fitted to.
std::vector< Int_t > fExcludedBins
Bins excluded from the fit.
TH1 * GetMCPrediction(Int_t parm) const
Return the adjusted MC template (Aji) for template (parm).
void CheckParNo(Int_t parm) const
Function for internal use, checking parameter validity An invalid parameter results in an error.
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:109
virtual void SetDirectory(TDirectory *dir)
By default, when a histogram is created, it is added to the list of histogram objects in the current ...
Definition TH1.cxx:8978
static TClass * Class()
virtual Int_t GetNbinsY() const
Definition TH1.h:543
virtual Int_t GetNbinsZ() const
Definition TH1.h:544
virtual Int_t GetDimension() const
Definition TH1.h:528
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition TH1.cxx:7151
virtual Int_t GetBin(Int_t binx, Int_t biny=0, Int_t binz=0) const
Return Global bin number corresponding to binx,y,z.
Definition TH1.cxx:4974
virtual Int_t GetNbinsX() const
Definition TH1.h:542
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition TH1.cxx:9260
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5076
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
Definition TH1.cxx:2734
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:173
const char * GetName() const override
Returns name of object.
Definition TNamed.h:49
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition TNamed.cxx:149
An array of TObjects.
Definition TObjArray.h:31
virtual void Expand(Int_t newSize)
Expand or shrink the array to newSize elements.
void AddAt(TObject *obj, Int_t idx) override
Add object at position ids.
void Delete(Option_t *option="") override
Remove all objects from the array AND delete all heap based objects.
TObject * At(Int_t idx) const override
Definition TObjArray.h:164
TObject * RemoveAt(Int_t idx) override
Remove object at index idx.
void Add(TObject *obj) override
Definition TObjArray.h:68
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition TObject.cxx:242
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:1057
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1071
Basic string class.
Definition TString.h:138
void ToUpper()
Change string to upper case.
Definition TString.cxx:1202
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2384
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:640
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
Int_t Nint(T x)
Round to nearest integer. Rounds half integers to the nearest even integer.
Definition TMath.h:704
Double_t Prob(Double_t chi2, Int_t ndf)
Computation of the probability for a certain Chi-squared (chi2) and number of degrees of freedom (ndf...
Definition TMath.cxx:637
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Definition TMath.h:767
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124