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