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