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
163
164
165////////////////////////////////////////////////////////////////////////////////
166/// TFractionFitter default constructor.
167
169fFitDone(kFALSE),
170fLowLimitX(0), fHighLimitX(0),
171fLowLimitY(0), fHighLimitY(0),
172fLowLimitZ(0), fHighLimitZ(0),
173fData(nullptr), fIntegralData(0),
174fPlot(nullptr)
175{
176 fFractionFitter = nullptr;
177 fIntegralMCs = nullptr;
178 fFractions = nullptr;
179
180 fNpfits = 0;
181 fNDF = 0;
182 fChisquare = 0;
183 fNpar = 0;
184}
185
186////////////////////////////////////////////////////////////////////////////////
187/// TFractionFitter constructor. Does a complete initialisation (including
188/// consistency checks, default fit range as the whole histogram but without
189/// under- and overflows, and declaration of the fit parameters). Note that
190/// the histograms are not copied, only references are used.
191/// \param[in] data histogram to be fitted
192/// \param[in] MCs array of TH1* corresponding template distributions
193/// \param[in] option can be used to control the print level of the minimization algorithm
194/// - option = "Q" : quite - no message is printed
195/// - option = "V" : verbose - max print out
196/// - option = "" : default: print initial fraction values and result
197
199fFitDone(kFALSE), fChisquare(0), fPlot(nullptr) {
200 fData = data;
201 // Default: include all of the histogram (but without under- and overflows)
202 fLowLimitX = 1;
204 if (fData->GetDimension() > 1) {
205 fLowLimitY = 1;
207 if (fData->GetDimension() > 2) {
208 fLowLimitZ = 1;
210 }
211 }
212 fNpar = MCs->GetEntries();
213 Int_t par;
214 for (par = 0; par < fNpar; ++par) {
215 fMCs.Add(MCs->At(par));
216 // Histogram containing template prediction
217 TString s = TString::Format("Prediction for MC sample %i",par);
218 TH1* pred = (TH1*) ((TH1*)MCs->At(par))->Clone(s);
219 // TFractionFitter manages these histograms
220 pred->SetDirectory(nullptr);
221 pred->SetTitle(s);
222 fAji.Add(pred);
223 }
226
229
231
232 // set print level
233 TString opt(option);
234 opt.ToUpper();
235 if (opt.Contains("Q") ) {
236 fFractionFitter->Config().MinimizerOptions().SetPrintLevel(0);
237 }
238 else if (opt.Contains("V") ) {
239 fFractionFitter->Config().MinimizerOptions().SetPrintLevel(2);
240 }
241 else
242 fFractionFitter->Config().MinimizerOptions().SetPrintLevel(1);
243
245 Double_t defaultStep = 0.01;
246 // set the parameters
247 std::vector<ROOT::Fit::ParameterSettings> & parameters = fFractionFitter->Config().ParamsSettings();
248 parameters.reserve(fNpar);
249 for (par = 0; par < fNpar; ++par) {
250 TString name("frac"); name += par;
251 parameters.push_back(ROOT::Fit::ParameterSettings(name.Data(), defaultFraction, defaultStep) );
252 }
253
254 if (fFractionFitter->Config().MinimizerOptions().ErrorDef() == 1.0 )
255 fFractionFitter->Config().MinimizerOptions().SetErrorDef(0.5);
256
257}
258
259////////////////////////////////////////////////////////////////////////////////
260/// TFractionFitter default destructor
261
264 delete[] fIntegralMCs;
265 delete[] fFractions;
266 if (fPlot) delete fPlot;
267 fAji.Delete();
268}
269
270////////////////////////////////////////////////////////////////////////////////
271/// Change the histogram to be fitted to. Notes:
272/// - Parameter constraints and settings are retained from a possible previous fit.
273/// - Modifying the dimension or number of bins results in an error (in this case
274/// rather instantiate a new TFractionFitter object)
275
281
282////////////////////////////////////////////////////////////////////////////////
283/// Change the histogram for template number `<parm>`. Notes:
284/// - Parameter constraints and settings are retained from a possible previous fit.
285/// - Modifying the dimension or number of bins results in an error (in this case
286/// rather instantiate a new TFractionFitter object)
287
295
296////////////////////////////////////////////////////////////////////////////////
297/// Set bin by bin weights for template number `<parm>` (the parameter numbering
298/// follows that of the input template vector).
299/// Weights can be "unset" by passing a null pointer.
300/// Consistency of the weights histogram with the data histogram is checked at
301/// this point, and an error in case of problems.
302
305 if (fWeights[parm]) {
307 }
308 if (weight) {
309 if (weight->GetNbinsX() != fData->GetNbinsX() ||
310 (fData->GetDimension() > 1 && weight->GetNbinsY() != fData->GetNbinsY()) ||
311 (fData->GetDimension() > 2 && weight->GetNbinsZ() != fData->GetNbinsZ())) {
312 Error("SetWeight","Inconsistent weights histogram for source %d", parm);
313 return;
314 }
315 TString ts = "weight hist: "; ts += weight->GetName();
316 fWeights.AddAt(weight,parm);
317 }
318}
319
320////////////////////////////////////////////////////////////////////////////////
321/// Give direct access to the underlying fitter class. This can be
322/// used e.g. to modify parameter values or step sizes.
323
327
328////////////////////////////////////////////////////////////////////////////////
329/// Function for internal use, checking parameter validity
330/// An invalid parameter results in an error.
331
334 Error("CheckParNo","Invalid parameter number %d",parm);
335 }
336}
337
338////////////////////////////////////////////////////////////////////////////////
339/// Set the X range of the histogram to be used in the fit.
340/// Use ReleaseRangeX() to go back to fitting the full histogram.
341/// The consistency check ensures that no empty fit range occurs (and also
342/// recomputes the bin content integrals).
343/// \param[in] low lower X bin number
344/// \param[in] high upper X bin number
345
347 fLowLimitX = (low > 0 ) ? low : 1;
348 fHighLimitX = ( high > 0 && high <= fData->GetNbinsX()) ? high : fData->GetNbinsX();
350}
351
352////////////////////////////////////////////////////////////////////////////////
353/// Release restrictions on the X range of the histogram to be used in the fit.
354
360
361////////////////////////////////////////////////////////////////////////////////
362/// Set the Y range of the histogram to be used in the fit (2D or 3D histograms only).
363/// Use ReleaseRangeY() to go back to fitting the full histogram.
364/// The consistency check ensures that no empty fit range occurs (and also
365/// recomputes the bin content integrals).
366/// \param[in] low lower X bin number
367/// \param[in] high upper X bin number
368
370 if (fData->GetDimension() < 2) {
371 Error("SetRangeY","Y range cannot be set for 1D histogram");
372 return;
373 }
374
375 fLowLimitY = (low > 0) ? low : 1;
376 fHighLimitY = (high > 0 && high <= fData->GetNbinsY()) ? high : fData->GetNbinsY();
378}
379
380////////////////////////////////////////////////////////////////////////////////
381/// Release restrictions on the Y range of the histogram to be used in the fit.
382
388
389
390////////////////////////////////////////////////////////////////////////////////
391/// Set the Z range of the histogram to be used in the fit (3D histograms only).
392/// Use ReleaseRangeY() to go back to fitting the full histogram.
393/// The consistency check ensures that no empty fit range occurs (and also
394/// recomputes the bin content integrals).
395/// \param[in] low lower X bin number
396/// \param[in] high upper X bin number
397
399 if (fData->GetDimension() < 3) {
400 Error("SetRangeZ","Z range cannot be set for 1D or 2D histogram");
401 return;
402 }
403
404
405 fLowLimitZ = (low > 0) ? low : 1;
406 fHighLimitZ = (high > 0 && high <= fData->GetNbinsZ()) ? high : fData->GetNbinsZ();
408}
409
410////////////////////////////////////////////////////////////////////////////////
411/// Release restrictions on the Z range of the histogram to be used in the fit.
412
418
419////////////////////////////////////////////////////////////////////////////////
420/// Exclude the given bin from the fit. The bin numbering to be used is that
421/// of TH1::GetBin().
422
424 int excluded = fExcludedBins.size();
425 for (int b = 0; b < excluded; ++b) {
426 if (fExcludedBins[b] == bin) {
427 Error("ExcludeBin", "bin %d already excluded", bin);
428 return;
429 }
430 }
431 fExcludedBins.push_back(bin);
432 // This call serves to properly (re)determine the number of degrees of freeom
434}
435
436////////////////////////////////////////////////////////////////////////////////
437/// Include the given bin in the fit, if it was excluded before using ExcludeBin().
438/// The bin numbering to be used is that of TH1::GetBin().
439
441 for (std::vector<Int_t>::iterator it = fExcludedBins.begin();
442 it != fExcludedBins.end(); ++it) {
443 if (*it == bin) {
444 fExcludedBins.erase(it);
445 // This call serves to properly (re)determine the number of degrees of freeom
447 return;
448 }
449 }
450 Error("IncludeBin", "bin %d was not excluded", bin);
451}
452
453////////////////////////////////////////////////////////////////////////////////
454/// Function for internal use, checking whether the given bin is
455/// excluded from the fit or not.
456
458 for (unsigned int b = 0; b < fExcludedBins.size(); ++b)
459 if (fExcludedBins[b] == bin) return true;
460 return false;
461}
462
463////////////////////////////////////////////////////////////////////////////////
464/// Constrain the values of parameter number `<parm>` (the parameter numbering
465/// follows that of the input template vector).
466/// Use UnConstrain() to remove this constraint.
467
470 assert( parm >= 0 && parm < (int) fFractionFitter->Config().ParamsSettings().size() );
471 fFractionFitter->Config().ParSettings(parm).SetLimits(low,high);
472}
473
474////////////////////////////////////////////////////////////////////////////////
475/// Remove the constraints on the possible values of parameter `<parm>`.
476
479 fFractionFitter->Config().ParSettings(parm).RemoveLimits();
480}
481
482////////////////////////////////////////////////////////////////////////////////
483/// Function used internally to check the consistency between the
484/// various histograms. Checks are performed on nonexistent or empty
485/// histograms, the precise histogram class, and the number of bins.
486/// In addition, integrals over the "allowed" bin ranges are computed.
487/// Any inconsistency results in a error.
488
490 if (! fData) {
491 Error("CheckConsistency","Nonexistent data histogram");
492 return;
493 }
495 Int_t x,y,z,par;
497 fIntegralData = 0;
498 fNpfits = 0;
499 for (z = minZ; z <= maxZ; ++z) {
500 for (y = minY; y <= maxY; ++y) {
501 for (x = minX; x <= maxX; ++x) {
502 if (IsExcluded(fData->GetBin(x, y, z))) continue;
503 fNpfits++;
505 }
506 }
507 }
508 if (fIntegralData <= 0) {
509 Error("CheckConsistency","Empty data histogram");
510 return;
511 }
512 TClass* cl = fData->Class();
513
514 fNDF = fNpfits - fNpar;
515
516 if (fNpar < 2) {
517 Error("CheckConsistency","Need at least two MC histograms");
518 return;
519 }
520
521 for (par = 0; par < fNpar; ++par) {
522 TH1 *h = (TH1*)fMCs.At(par);
523 if (! h) {
524 Error("CheckConsistency","Nonexistent MC histogram for source #%d",par);
525 return;
526 }
527 if ((! h->Class()->InheritsFrom(cl)) || h->GetNbinsX() != fData->GetNbinsX() ||
528 (fData->GetDimension() > 1 && h->GetNbinsY() != fData->GetNbinsY()) ||
529 (fData->GetDimension() > 2 && h->GetNbinsZ() != fData->GetNbinsZ())) {
530 Error("CheckConsistency","Histogram inconsistency for source #%d",par);
531 return;
532 }
533 fIntegralMCs[par] = 0;
534 for (z = minZ; z <= maxZ; ++z) {
535 for (y = minY; y <= maxY; ++y) {
536 for (x = minX; x <= maxX; ++x) {
537 Int_t bin = fData->GetBin(x, y, z);
538 if (IsExcluded(bin)) continue;
539 Double_t MCEvents = h->GetBinContent(bin);
540 if (MCEvents < 0) {
541 Error("CheckConsistency", "Number of MC events (bin = %d, par = %d) cannot be negative: "
542 " their distribution is binomial (see paper)", bin, par);
543 }
544 fIntegralMCs[par] += MCEvents;
545 }
546 }
547 }
548 if (fIntegralMCs[par] <= 0) {
549 Error("CheckConsistency","Empty MC histogram #%d",par);
550 }
551 }
552}
553
554////////////////////////////////////////////////////////////////////////////////
555/// Perform the fit with the default UP value.
556/// The value returned is the minimisation status.
557
559
560 // remove any existing output histogram
561 if (fPlot) {
562 delete fPlot; fPlot = nullptr;
563 }
564
565 // Make sure the correct likelihood computation is used
568
569 // fit
570 Bool_t status = fFractionFitter->FitFCN();
571 if (!status) Warning("Fit","Abnormal termination of minimization.");
572
573 fFitDone = kTRUE;
574
575 // determine goodness of fit
577
578 // create a new result class
580 TString name = TString::Format("TFractionFitter_result_of_%s",fData->GetName() );
581 fr->SetName(name); fr->SetTitle(name);
582 return TFitResultPtr(fr);
583}
584
585////////////////////////////////////////////////////////////////////////////////
586/// Set UP to the given value (see class TMinuit), and perform a MINOS minimisation.
587
589 if (! fFitDone) {
590 Error("ErrorAnalysis","Fit not yet performed");
591 return;
592 }
593
594
595 Double_t up = UP > 0 ? UP : 0.5;
596 fFractionFitter->Config().MinimizerOptions().SetErrorDef(up);
598 if (!status) {
599 Error("ErrorAnalysis","Error return from MINOS: %d",fFractionFitter->Result().Status());
600 }
601}
602
603////////////////////////////////////////////////////////////////////////////////
604/// Obtain the fit result for parameter `<parm>` (the parameter numbering
605/// follows that of the input template vector).
606
609 if (! fFitDone) {
610 Error("GetResult","Fit not yet performed");
611 return;
612 }
613 value = fFractionFitter->Result().Parameter(parm);
614 error = fFractionFitter->Result().Error(parm);
615}
616
617////////////////////////////////////////////////////////////////////////////////
618/// Return the "template prediction" corresponding to the fit result (this is not
619/// the same as the weighted sum of template distributions, as template statistical
620/// uncertainties are taken into account).
621/// Note that the name of this histogram will simply be the same as that of the
622/// "data" histogram, prefixed with the string "Fraction fit to hist: ".
623/// Note also that the histogram is managed by the TFractionFitter class, so the returned pointer will be invalid if
624/// the class is deleted
625
627 if (! fFitDone) {
628 Error("GetPlot","Fit not yet performed");
629 return nullptr;
630 }
631 if (! fPlot) {
632 Double_t f = 0;
633 const Double_t * par = fFractionFitter->Result().GetParams();
634 assert(par);
635 ComputeFCN(f, par, 3);
636 }
637 return fPlot;
638}
639
640////////////////////////////////////////////////////////////////////////////////
641/// Used internally to obtain the bin ranges according to the dimensionality of
642/// the histogram and the limits set by hand.
643
645 Int_t& minZ, Int_t& maxZ) const {
646 if (fData->GetDimension() < 2) {
647 minY = maxY = minZ = maxZ = 0;
650 } else if (fData->GetDimension() < 3) {
651 minZ = maxZ = 0;
656 } else {
663 }
664}
665
666////////////////////////////////////////////////////////////////////////////////
667/// Used internally to compute the likelihood value.
668
670{
671 // normalise the fit parameters
672 Int_t bin, mc;
674 Int_t x,y,z;
676 for (mc = 0; mc < fNpar; ++mc) {
678 TH1 *h = (TH1*)fMCs[mc];
679 TH1 *hw = (TH1*)fWeights[mc];
680 if (hw) {
681 tot = 0;
682 for (z = minZ; z <= maxZ; ++z) {
683 for (y = minY; y <= maxY; ++y) {
684 for (x = minX; x <= maxX; ++x) {
685 if (IsExcluded(fData->GetBin(x, y, z))) continue;
686 Double_t weight = hw->GetBinContent(x, y, z);
687 if (weight <= 0) {
688 Error("ComputeFCN","Invalid weight encountered for MC source %d",mc);
689 return;
690 }
691 tot += weight * h->GetBinContent(x, y, z);
692 }
693 }
694 }
695 } else tot = fIntegralMCs[mc];
697 }
698
699 if (flag == 3) {
700 TString ts = "Fraction fit to hist: "; ts += fData->GetName();
701 fPlot = (TH1*) fData->Clone(ts.Data());
702 // plot histogram is managed by TFractionFitter
703 fPlot->SetDirectory(nullptr);
704 fPlot->Reset();
705 }
706 // likelihood computation
707 Double_t result = 0;
708 for (z = minZ; z <= maxZ; ++z) {
709 for (y = minY; y <= maxY; ++y) {
710 for (x = minX; x <= maxX; ++x) {
711 bin = fData->GetBin(x, y, z);
712 if (IsExcluded(bin)) continue;
713
714 // Solve for the "predictions"
715 int k0 = 0;
716 Double_t ti = 0.0; Double_t aki = 0.0;
717 FindPrediction(bin, ti, k0, aki);
718
720 for (mc = 0; mc < fNpar; ++mc) {
721 TH1 *h = (TH1*)fMCs[mc];
722 TH1 *hw = (TH1*)fWeights[mc];
724 Double_t binContent = h->GetBinContent(bin);
725 Double_t weight = hw ? hw->GetBinContent(bin) : 1;
726 if (k0 >= 0 && fFractions[mc] == fFractions[k0]) {
728 } else {
729 binPrediction = binContent > 0 ? binContent / (1+weight*fFractions[mc]*ti) : 0;
730 }
731
734 if (binContent > 0 && binPrediction > 0)
736
737 if (flag == 3) {
738 ((TH1*)fAji.At(mc))->SetBinContent(bin, binPrediction);
739 }
740 }
741
742 if (flag == 3) {
744 }
745
747 Double_t found = fData->GetBinContent(bin);
748 if (found > 0 && prediction > 0)
749 result += found*TMath::Log(prediction);
750 }
751 }
752 }
753
754 f = -result;
755}
756
757////////////////////////////////////////////////////////////////////////////////
758/// Function used internally to obtain the template prediction in the individual bins
759/// 'bin' <=> 'i' (paper)
760/// 'par' <=> 'j' (paper)
761
763 std::vector<Double_t> wgtFrac(fNpar); // weighted fractions (strengths of the sources)
764 std::vector<Double_t> a_ji(fNpar); // number of observed MC events for bin 'i' and par (source) 'j'
765 Double_t d_i = fData->GetBinContent(bin); // number of events in the real data for bin 'i'
766
767 // Cache the weighted fractions and the number of observed MC events
768 // Sanity check: none of the fractions should be == 0
769 for (Int_t par = 0; par < fNpar; ++par) {
770 a_ji[par] = ((TH1*)fMCs.At(par))->GetBinContent(bin);
771 TH1* hw = (TH1*)fWeights.At(par);
772 wgtFrac[par] = hw ? hw->GetBinContent(bin) * fFractions[par] : fFractions[par];
773 if (wgtFrac[par] == 0) {
774 Error("FindPrediction", "Fraction[%d] = 0!", par);
775 return;
776 }
777 }
778
779 // Case data = 0
780 if (TMath::Nint(d_i) == 0) {
781 t_i = 1;
782 k_0 = -1;
783 A_ki = 0;
784 return;
785 }
786
787 // Case one or more of the MC bin contents == 0 -> find largest fraction
788 // k_0 stores the source index of the largest fraction
789 k_0 = 0;
791 for (Int_t par = 1; par < fNpar; ++par) {
792 if (wgtFrac[par] > maxWgtFrac) {
793 k_0 = par;
794 maxWgtFrac = wgtFrac[par];
795 }
796 }
797 Double_t t_min = -1 / maxWgtFrac; // t_i cannot be smaller than this value (see paper, par 5)
798
799 // Determine if there are more sources which have the same maximum contribution (fraction)
801 for (Int_t par = 0; par < fNpar; ++par) {
802 if (par == k_0) continue;
803 if (wgtFrac[par] == maxWgtFrac) {
804 nMax++;
805 contentsMax += a_ji[par];
806 }
807 }
808
809 // special action if there is a zero in the number of entries for the MC source with
810 // the largest strength (fraction) -> See Paper, Paragraph 5
811 if (contentsMax == 0) {
812 A_ki = d_i / (1.0 + maxWgtFrac);
813 for (Int_t par = 0; par < fNpar; ++par) {
814 if (par == k_0 || wgtFrac[par] == maxWgtFrac) continue;
815 A_ki -= a_ji[par] * wgtFrac[par] / (maxWgtFrac - wgtFrac[par]);
816 }
817 if (A_ki > 0) {
818 A_ki /= nMax;
819 t_i = t_min;
820 return;
821 }
822 }
823 k_0 = -1;
824
825 // Case of nonzero histogram contents: solve for t_i using Newton's method
826 // The equation that needs to be solved:
827 // func(t_i) = \sum\limits_j{\frac{ p_j a_{ji} }{1 + p_j t_i}} - \frac{d_i}{1 - t_i} = 0
828 t_i = 0; Double_t step = 0.2;
829 Int_t maxIter = 100000; // maximum number of iterations
830 for(Int_t i = 0; i < maxIter; ++i) {
831 if (t_i >= 1 || t_i < t_min) {
832 step /= 10;
833 t_i = 0;
834 }
835 Double_t func = - d_i / (1.0 - t_i);
836 Double_t deriv = func / (1.0 - t_i);
837 for (Int_t par = 0; par < fNpar; ++par) {
838 Double_t r = 1.0 / (t_i + 1.0 / wgtFrac[par]);
839 func += a_ji[par] * r;
840 deriv -= a_ji[par] * r * r;
841 }
842 if (TMath::Abs(func) < 1e-12) return; // solution found
843 Double_t delta = - func / deriv; // update delta
844 if (TMath::Abs(delta) > step)
845 delta = (delta > 0) ? step : -step; // correct delta if it becomes too large
846 t_i += delta;
847 if (TMath::Abs(delta) < 1e-13) return; // solution found
848 } // the loop breaks when the solution is found, or when the maximum number of iterations is exhausted
849
850 Warning("FindPrediction", "Did not find solution for t_i in %d iterations", maxIter);
851
852 return;
853}
854
855#ifdef OLD
856////////////////////////////////////////////////////////////////////////////////
857/// Function called by the minimisation package. The actual functionality is passed
858/// on to the TFractionFitter::ComputeFCN member function.
859
861 TFractionFitter* fitter = dynamic_cast<TFractionFitter*>(fFractionFitter->GetObjectFit());
862 if (!fitter) {
863 Error("TFractionFitFCN","Invalid fit object encountered!");
864 return;
865 }
866 fitter->ComputeFCN(npar, gin, f, par, flag);
867}
868#endif
869
870////////////////////////////////////////////////////////////////////////////////
871/// Return the likelihood ratio Chi-squared (chi2) for the fit.
872/// The value is computed when the fit is executed successfully.
873/// Chi2 calculation is based on the "likelihood ratio" lambda,
874/// lambda = L(y;n) / L(m;n),
875/// where L(y;n) is the likelihood of the fit result `<y>` describing the data `<n>`
876/// and L(m;n) is the likelihood of an unknown "true" underlying distribution
877/// `<m>` describing the data `<n>`. Since `<m>` is unknown, the data distribution is
878/// used instead,
879/// lambda = L(y;n) / L(n;n).
880/// Note that this ratio is 1 if the fit is perfect. The chi2 value is then
881/// computed according to
882/// chi2 = -2*ln(lambda).
883/// This parameter can be shown to follow a Chi-square distribution. See for
884/// example S. Baker and R. Cousins, "Clarification of the use of chi-square
885/// and likelihood functions in fits to histograms", Nucl. Instr. Meth. A221,
886/// pp. 437-442 (1984)
887
892
893////////////////////////////////////////////////////////////////////////////////
894/// return the number of degrees of freedom in the fit
895/// the fNDF parameter has been previously computed during a fit.
896/// The number of degrees of freedom corresponds to the number of points
897/// used in the fit minus the number of templates.
898
900{
901 if (fNDF == 0) return fNpfits-fNpar;
902 return fNDF;
903}
904
905////////////////////////////////////////////////////////////////////////////////
906/// return the fit probability
907
909{
911 if (ndf <= 0) return 0;
912 return TMath::Prob(fChisquare,ndf);
913}
914
915////////////////////////////////////////////////////////////////////////////////
916/// Method used internally to compute the likelihood ratio chi2
917/// See the function GetChisquare() for details
918
920{
921 if ( !fFitDone ) {
922 Error("ComputeChisquareLambda","Fit not yet (successfully) performed");
923 fChisquare = 0;
924 return;
925 }
926
927 // fPlot must be initialized and filled. Leave this to the GetPlot() method.
928 if (! fPlot)
929 GetPlot();
930
933
934 Double_t logLyn = 0; // likelihood of prediction
935 Double_t logLmn = 0; // likelihood of data ("true" distribution)
936 for(Int_t x = minX; x <= maxX; x++) {
937 for(Int_t y = minY; y <= maxY; y++) {
938 for(Int_t z = minZ; z <= maxZ; z++) {
939 if (IsExcluded(fData->GetBin(x, y, z))) continue;
942 if(fi != 0) logLyn += di * TMath::Log(fi) - fi;
943 if(di != 0) logLmn += di * TMath::Log(di) - di;
944 for(Int_t j = 0; j < fNpar; j++) {
945 Double_t aji = ((TH1*)fMCs.At(j))->GetBinContent(x, y, z);
946 Double_t bji = ((TH1*)fAji.At(j))->GetBinContent(x, y, z);
947 if(bji != 0) logLyn += aji * TMath::Log(bji) - bji;
948 if(aji != 0) logLmn += aji * TMath::Log(aji) - aji;
949 }
950 }
951 }
952 }
953
954 fChisquare = -2*logLyn + 2*logLmn;
955
956 return;
957}
958
959////////////////////////////////////////////////////////////////////////////////
960/// Return the adjusted MC template (Aji) for template (parm).
961/// Note that the (Aji) times fractions only sum to the total prediction
962/// of the fit if all weights are 1.
963/// Note also that the histogram is managed by the TFractionFitter class, so the returned pointer will be invalid if
964/// the class is deleted
965
967{
969 if ( !fFitDone ) {
970 Error("GetMCPrediction","Fit not yet performed");
971 return nullptr;
972 }
973 return (TH1*) fAji.At(parm);
974}
#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
#define ClassImp(name)
Definition Rtypes.h:376
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
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:8966
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:7139
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:4964
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:9248
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5064
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
Definition TH1.cxx:2724
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:174
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:150
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:243
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:1058
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1072
Basic string class.
Definition TString.h:138
void ToUpper()
Change string to upper case.
Definition TString.cxx:1203
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:2385
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:699
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:762
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124