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