1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 /* $Id: AliUnfolding.cxx 31168 2009-02-23 15:18:45Z jgrosseo $ */
18 // This class allows 1-dimensional unfolding.
19 // Methods that are implemented are chi2 minimization and bayesian unfolding.
21 // Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
23 #include "AliUnfolding.h"
26 #include <TVirtualFitter.h>
31 #include "Riostream.h"
34 using namespace std; //required for resolving the 'cout' symbol
36 TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
37 TMatrixD* AliUnfolding::fgCorrelationMatrixSquared = 0;
38 TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
39 TVectorD* AliUnfolding::fgCurrentESDVector = 0;
40 TVectorD* AliUnfolding::fgEntropyAPriori = 0;
41 TVectorD* AliUnfolding::fgEfficiency = 0;
43 TAxis* AliUnfolding::fgUnfoldedAxis = 0;
44 TAxis* AliUnfolding::fgMeasuredAxis = 0;
46 TF1* AliUnfolding::fgFitFunction = 0;
48 AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid;
49 Int_t AliUnfolding::fgMaxInput = -1; // bins in measured histogram
50 Int_t AliUnfolding::fgMaxParams = -1; // bins in unfolded histogram = number of fit params
51 Float_t AliUnfolding::fgOverflowBinLimit = -1;
53 AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1;
54 Float_t AliUnfolding::fgRegularizationWeight = 10000;
55 Int_t AliUnfolding::fgSkipBinsBegin = 0;
56 Float_t AliUnfolding::fgMinuitStepSize = 0.1; // (usually not needed to be changed) step size in minimization
57 Float_t AliUnfolding::fgMinuitPrecision = 1e-6; // minuit precision
58 Int_t AliUnfolding::fgMinuitMaxIterations = 1e6; // minuit maximum number of iterations
59 Double_t AliUnfolding::fgMinuitStrategy = 1.; // minuit strategy
60 Bool_t AliUnfolding::fgMinimumInitialValue = kFALSE; // set all initial values at least to the smallest value among the initial values
61 Float_t AliUnfolding::fgMinimumInitialValueFix = -1;
62 Bool_t AliUnfolding::fgNormalizeInput = kFALSE; // normalize input spectrum
63 Float_t AliUnfolding::fgNotFoundEvents = 0;
64 Bool_t AliUnfolding::fgSkipBin0InChi2 = kFALSE;
66 Float_t AliUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
67 Int_t AliUnfolding::fgBayesianIterations = 10; // number of iterations in Bayesian method
69 Bool_t AliUnfolding::fgDebug = kFALSE;
71 Int_t AliUnfolding::fgCallCount = 0;
73 Int_t AliUnfolding::fgPowern = 5;
75 Double_t AliUnfolding::fChi2FromFit = 0.;
76 Double_t AliUnfolding::fPenaltyVal = 0.;
77 Double_t AliUnfolding::fAvgResidual = 0.;
79 Int_t AliUnfolding::fgPrintChi2Details = 0;
81 TCanvas *AliUnfolding::fgCanvas = 0;
82 TH1 *AliUnfolding::fghUnfolded = 0;
83 TH2 *AliUnfolding::fghCorrelation = 0;
84 TH1 *AliUnfolding::fghEfficiency = 0;
85 TH1 *AliUnfolding::fghMeasured = 0;
87 ClassImp(AliUnfolding)
89 //____________________________________________________________________
90 void AliUnfolding::SetUnfoldingMethod(MethodType methodType)
92 // set unfolding method
93 fgMethodType = methodType;
98 case kInvalid: name = "INVALID"; break;
99 case kChi2Minimization: name = "Chi2 Minimization"; break;
100 case kBayesian: name = "Bayesian unfolding"; break;
101 case kFunction: name = "Functional fit"; break;
103 Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name);
106 //____________________________________________________________________
107 void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit)
109 // enable the creation of a overflow bin that includes all statistics below the given limit
111 fgOverflowBinLimit = overflowBinLimit;
113 Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit);
116 //____________________________________________________________________
117 void AliUnfolding::SetSkipBinsBegin(Int_t nBins)
119 // set number of skipped bins in regularization
121 fgSkipBinsBegin = nBins;
123 Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin);
126 //____________________________________________________________________
127 void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded)
129 // set number of bins in the input (measured) distribution and in the unfolded distribution
130 fgMaxInput = nMeasured;
131 fgMaxParams = nUnfolded;
133 if (fgCorrelationMatrix)
135 delete fgCorrelationMatrix;
136 fgCorrelationMatrix = 0;
138 if (fgCorrelationMatrixSquared)
140 fgCorrelationMatrixSquared = 0;
141 delete fgCorrelationMatrixSquared;
143 if (fgCorrelationCovarianceMatrix)
145 delete fgCorrelationCovarianceMatrix;
146 fgCorrelationCovarianceMatrix = 0;
148 if (fgCurrentESDVector)
150 delete fgCurrentESDVector;
151 fgCurrentESDVector = 0;
153 if (fgEntropyAPriori)
155 delete fgEntropyAPriori;
156 fgEntropyAPriori = 0;
165 delete fgUnfoldedAxis;
170 delete fgMeasuredAxis;
174 Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
177 //____________________________________________________________________
178 void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight)
181 // sets the parameters for chi2 minimization
184 fgRegularizationType = type;
185 fgRegularizationWeight = weight;
187 Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight);
190 //____________________________________________________________________
191 void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
194 // sets the parameters for Bayesian unfolding
197 fgBayesianSmoothing = smoothing;
198 fgBayesianIterations = nIterations;
200 Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
203 //____________________________________________________________________
204 void AliUnfolding::SetFunction(TF1* function)
206 // set function for unfolding with a fit function
208 fgFitFunction = function;
210 Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar());
213 //____________________________________________________________________
214 Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
216 // unfolds with unfolding method fgMethodType
219 // correlation: response matrix as measured vs. generated
220 // efficiency: (optional) efficiency that is applied on the unfolded spectrum, i.e. it has to be in unfolded variables. If 0 no efficiency is applied.
221 // measured: the measured spectrum
222 // initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions.
223 // result: target for the unfolded result
224 // check: depends on the unfolding method, see comments in specific functions
226 // return code: see UnfoldWithMinuit/UnfoldWithBayesian/UnfoldWithFunction
228 if (fgMaxInput == -1)
230 Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
231 fgMaxInput = measured->GetNbinsX();
233 if (fgMaxParams == -1)
235 Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
236 fgMaxParams = measured->GetNbinsX();
239 if (fgOverflowBinLimit > 0)
240 CreateOverflowBin(correlation, measured);
242 switch (fgMethodType)
246 Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting...");
249 case kChi2Minimization:
250 return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check);
252 return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result);
254 return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result);
262 //____________________________________________________________________
263 void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency)
265 // fill static variables needed for minuit fit
267 if (!fgCorrelationMatrix)
268 fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
269 if (!fgCorrelationMatrixSquared)
270 fgCorrelationMatrixSquared = new TMatrixD(fgMaxInput, fgMaxParams);
271 if (!fgCorrelationCovarianceMatrix)
272 fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
273 if (!fgCurrentESDVector)
274 fgCurrentESDVector = new TVectorD(fgMaxInput);
275 if (!fgEntropyAPriori)
276 fgEntropyAPriori = new TVectorD(fgMaxParams);
278 fgEfficiency = new TVectorD(fgMaxParams);
280 delete fgUnfoldedAxis;
281 fgUnfoldedAxis = new TAxis(*(correlation->GetXaxis()));
283 delete fgMeasuredAxis;
284 fgMeasuredAxis = new TAxis(*(correlation->GetYaxis()));
286 fgCorrelationMatrix->Zero();
287 fgCorrelationCovarianceMatrix->Zero();
288 fgCurrentESDVector->Zero();
289 fgEntropyAPriori->Zero();
291 // normalize correction for given nPart
292 for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
294 Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
297 Float_t maxValue = 0;
299 for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
301 // find most probably value
302 if (maxValue < correlation->GetBinContent(i, j))
304 maxValue = correlation->GetBinContent(i, j);
309 correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i));
310 correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
312 if (i <= fgMaxParams && j <= fgMaxInput)
314 (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
315 (*fgCorrelationMatrixSquared)(j-1, i-1) = correlation->GetBinContent(i, j) * correlation->GetBinContent(i, j);
319 //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
323 Float_t smallestError = 1;
324 if (fgNormalizeInput)
326 Float_t sumMeasured = measured->Integral();
327 measured->Scale(1.0 / sumMeasured);
328 smallestError /= sumMeasured;
331 for (Int_t i=0; i<fgMaxInput; ++i)
333 (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
334 if (measured->GetBinError(i+1) > 0)
336 (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
338 else // in this case put error of 1, otherwise 0 bins are not added to the chi2...
339 (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / smallestError / smallestError;
341 if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
342 (*fgCorrelationCovarianceMatrix)(i, i) = 0;
343 //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
346 // efficiency is expected to match bin width of result
347 for (Int_t i=0; i<fgMaxParams; ++i)
349 (*fgEfficiency)(i) = efficiency->GetBinContent(i+1);
352 if (correlation->GetNbinsX() != fgMaxParams || correlation->GetNbinsY() != fgMaxInput)
353 cout << "Response histo has incorrect dimensions; expect (" << fgMaxParams << ", " << fgMaxInput << "), got (" << correlation->GetNbinsX() << ", " << correlation->GetNbinsY() << ")" << endl;
357 //____________________________________________________________________
358 Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
361 // implementation of unfolding (internal function)
363 // unfolds <measured> using response from <correlation> and effiency <efficiency>
364 // output is in <result>
365 // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
366 // negative values in initialConditions mean that the given parameter is fixed to the absolute of the value
367 // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
369 // returns minuit status (0 = success), (-1 when check was set)
372 SetStaticVariables(correlation, measured, efficiency);
374 // Initialize TMinuit via generic fitter interface
375 Int_t params = fgMaxParams;
376 if (fgNotFoundEvents > 0)
379 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, params);
380 Double_t arglist[100];
381 // minuit->SetDefaultFitter("Minuit2");
383 // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
385 minuit->ExecuteCommand("SET PRINT", arglist, 1);
387 // however, enable warnings
388 //minuit->ExecuteCommand("SET WAR", arglist, 0);
390 // set minimization function
391 minuit->SetFCN(Chi2Function);
394 minuit->SetPrecision(fgMinuitPrecision);
396 minuit->SetMaxIterations(fgMinuitMaxIterations);
398 minuit->ExecuteCommand("SET STRATEGY",&fgMinuitStrategy,1);
400 for (Int_t i=0; i<fgMaxParams; i++)
401 (*fgEntropyAPriori)[i] = 1;
403 // set initial conditions as a-priori distribution for MRX regularization
405 for (Int_t i=0; i<fgMaxParams; i++)
406 if (initialConditions && initialConditions->GetBinContent(i+1) > 0)
407 (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
410 if (!initialConditions) {
411 initialConditions = measured;
413 Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions...");
414 //new TCanvas; initialConditions->DrawCopy();
415 if (fgNormalizeInput)
416 initialConditions->Scale(1.0 / initialConditions->Integral());
419 // extract minimum value from initial conditions (if we set a value to 0 it will stay 0)
420 Float_t minValue = 1e35;
421 if (fgMinimumInitialValueFix < 0)
423 for (Int_t i=0; i<fgMaxParams; ++i)
425 Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
426 if (initialConditions->GetBinContent(bin) > 0)
427 minValue = TMath::Min(minValue, (Float_t) initialConditions->GetBinContent(bin));
431 minValue = fgMinimumInitialValueFix;
433 Double_t* results = new Double_t[fgMaxParams+1];
434 for (Int_t i=0; i<fgMaxParams; ++i)
436 Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
437 results[i] = initialConditions->GetBinContent(bin);
443 results[i] = -results[i];
446 if (!fix && fgMinimumInitialValue && results[i] < minValue)
447 results[i] = minValue;
449 // minuit sees squared values to prevent it from going negative...
450 results[i] = TMath::Sqrt(results[i]);
452 minuit->SetParameter(i, Form("param%d", i), results[i], (fix) ? 0 : fgMinuitStepSize, 0, 0);
454 if (fgNotFoundEvents > 0)
456 results[fgMaxParams] = efficiency->GetBinContent(1);
457 minuit->SetParameter(fgMaxParams, "vtx0", results[fgMaxParams], fgMinuitStepSize / 100, 0.01, 0.80);
462 Chi2Function(dummy, 0, chi2, results, 0);
463 printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2);
472 // first param is number of iterations, second is precision....
473 arglist[0] = (float)fgMinuitMaxIterations;
474 // arglist[1] = 1e-5;
475 // minuit->ExecuteCommand("SET PRINT", arglist, 3);
476 // minuit->ExecuteCommand("SCAN", arglist, 0);
477 Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
478 Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status);
479 //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
480 //minuit->ExecuteCommand("MINOS", arglist, 0);
482 if (fgNotFoundEvents > 0)
484 results[fgMaxParams] = minuit->GetParameter(fgMaxParams);
485 Printf("Efficiency for bin 0 changed from %f to %f", efficiency->GetBinContent(1), results[fgMaxParams]);
486 efficiency->SetBinContent(1, results[fgMaxParams]);
489 for (Int_t i=0; i<fgMaxParams; ++i)
491 results[i] = minuit->GetParameter(i);
492 Double_t value = results[i] * results[i];
493 // error is : 2 * (relError on results[i]) * (value) = 2 * (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
495 if (TMath::IsNaN(minuit->GetParError(i)))
496 Printf("WARNING: Parameter %d error is nan", i);
498 error = 2 * minuit->GetParError(i) * results[i];
502 //printf("value before efficiency correction: %f\n",value);
503 if (efficiency->GetBinContent(i+1) > 0)
505 value /= efficiency->GetBinContent(i+1);
506 error /= efficiency->GetBinContent(i+1);
514 //printf("value after efficiency correction: %f +/- %f\n",value,error);
515 result->SetBinContent(i+1, value);
516 result->SetBinError(i+1, error);
519 Int_t tmpCallCount = fgCallCount;
520 fgCallCount = 0; // needs to be 0 so that the Chi2Function prints its output
521 Chi2Function(dummy, 0, chi2, results, 0);
523 Printf("AliUnfolding::UnfoldWithMinuit: iterations %d. Chi2 of final parameters is = %f", tmpCallCount, chi2);
530 //____________________________________________________________________
531 Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult)
534 // unfolds a spectrum using the Bayesian method
537 if (measured->Integral() <= 0)
539 Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
543 const Int_t kStartBin = 0;
545 Int_t kMaxM = fgMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
546 Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
548 // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
549 const Double_t kConvergenceLimit = kMaxT * 1e-6;
551 // store information in arrays, to increase processing speed (~ factor 5)
552 Double_t* measuredCopy = new Double_t[kMaxM];
553 Double_t* measuredError = new Double_t[kMaxM];
554 Double_t* prior = new Double_t[kMaxT];
555 Double_t* result = new Double_t[kMaxT];
556 Double_t* efficiency = new Double_t[kMaxT];
557 Double_t* binWidths = new Double_t[kMaxT];
559 Double_t** response = new Double_t*[kMaxT];
560 Double_t** inverseResponse = new Double_t*[kMaxT];
561 for (Int_t i=0; i<kMaxT; i++)
563 response[i] = new Double_t[kMaxM];
564 inverseResponse[i] = new Double_t[kMaxM];
568 Float_t measuredIntegral = measured->Integral();
569 for (Int_t m=0; m<kMaxM; m++)
571 measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
572 measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
574 for (Int_t t=0; t<kMaxT; t++)
576 response[t][m] = correlation->GetBinContent(t+1, m+1);
577 inverseResponse[t][m] = 0;
581 for (Int_t t=0; t<kMaxT; t++)
585 efficiency[t] = aEfficiency->GetBinContent(t+1);
590 prior[t] = measuredCopy[t];
592 binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1);
595 // pick prior distribution
596 if (initialConditions)
598 printf("Using different starting conditions...\n");
600 Float_t inputDistIntegral = initialConditions->Integral();
601 for (Int_t i=0; i<kMaxT; i++)
602 prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
605 //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
609 for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations < 0; i++)
612 Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
614 // calculate IR from Bayes theorem
615 // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
617 Double_t chi2Measured = 0;
618 for (Int_t m=0; m<kMaxM; m++)
621 for (Int_t t = kStartBin; t<kMaxT; t++)
622 norm += response[t][m] * prior[t];
624 // calc. chi2: (measured - response * prior) / error
625 if (measuredError[m] > 0)
627 Double_t value = (measuredCopy[m] - norm) / measuredError[m];
628 chi2Measured += value * value;
633 for (Int_t t = kStartBin; t<kMaxT; t++)
634 inverseResponse[t][m] = response[t][m] * prior[t] / norm;
638 for (Int_t t = kStartBin; t<kMaxT; t++)
639 inverseResponse[t][m] = 0;
642 //Printf("chi2Measured of the last prior is %e", chi2Measured);
644 for (Int_t t = kStartBin; t<kMaxT; t++)
647 for (Int_t m=0; m<kMaxM; m++)
648 value += inverseResponse[t][m] * measuredCopy[m];
650 if (efficiency[t] > 0)
651 result[t] = value / efficiency[t];
657 // draw intermediate result
658 for (Int_t t=0; t<kMaxT; t++)
660 aResult->SetBinContent(t+1, result[t]);
662 aResult->SetMarkerStyle(24+i);
663 aResult->SetMarkerColor(2);
664 aResult->DrawCopy((i == 0) ? "P" : "PSAME");
667 Double_t chi2LastIter = 0;
668 // regularization (simple smoothing)
669 for (Int_t t=kStartBin; t<kMaxT; t++)
671 Float_t newValue = 0;
673 // 0 bin excluded from smoothing
674 if (t > kStartBin+2 && t<kMaxT-1)
676 Float_t average = (result[t-1] / binWidths[t-1] + result[t] / binWidths[t] + result[t+1] / binWidths[t+1]) / 3 * binWidths[t];
678 // weight the average with the regularization parameter
679 newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
682 newValue = result[t];
684 // calculate chi2 (change from last iteration)
687 Double_t diff = (prior[t] - newValue) / prior[t];
688 chi2LastIter += diff * diff;
693 //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
694 //convergence->Fill(i+1, chi2LastIter);
696 if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit)
698 Printf("AliUnfolding::UnfoldWithBayesian: Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured);
701 } // end of iterations
703 //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
704 //delete convergence;
707 if (!fgNormalizeInput)
708 factor = measuredIntegral;
709 for (Int_t t=0; t<kMaxT; t++)
710 aResult->SetBinContent(t+1, result[t] * factor);
712 delete[] measuredCopy;
713 delete[] measuredError;
719 for (Int_t i=0; i<kMaxT; i++)
721 delete[] response[i];
722 delete[] inverseResponse[i];
725 delete[] inverseResponse;
730 // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
732 /*printf("Calculating covariance matrix. This may take some time...\n");
734 // check if this is the right one...
735 TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
737 Int_t xBins = hInverseResponseBayes->GetNbinsX();
738 Int_t yBins = hInverseResponseBayes->GetNbinsY();
740 // calculate "unfolding matrix" Mij
741 Float_t matrixM[251][251];
742 for (Int_t i=1; i<=xBins; i++)
744 for (Int_t j=1; j<=yBins; j++)
746 if (fCurrentEfficiency->GetBinContent(i) > 0)
747 matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
749 matrixM[i-1][j-1] = 0;
753 Float_t* vectorn = new Float_t[yBins];
754 for (Int_t j=1; j<=yBins; j++)
755 vectorn[j-1] = fCurrentESD->GetBinContent(j);
757 // first part of covariance matrix, depends on input distribution n(E)
758 Float_t cov1[251][251];
760 Float_t nEvents = fCurrentESD->Integral(); // N
765 for (Int_t k=0; k<xBins; k++)
767 printf("In Cov1: %d\n", k);
768 for (Int_t l=0; l<yBins; l++)
772 // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
773 for (Int_t j=0; j<yBins; j++)
774 cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
775 * (1.0 - vectorn[j] / nEvents);
777 // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
778 for (Int_t i=0; i<yBins; i++)
779 for (Int_t j=0; j<yBins; j++)
783 cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
784 * vectorn[j] / nEvents;
789 printf("Cov1 finished\n");
791 TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
794 for (Int_t i=1; i<=xBins; i++)
795 for (Int_t j=1; j<=yBins; j++)
796 cov->SetBinContent(i, j, cov1[i-1][j-1]);
801 // second part of covariance matrix, depends on response matrix
802 Float_t cov2[251][251];
804 // Cov[P(Er|Cu), P(Es|Cu)] term
805 Float_t covTerm[100][100][100];
806 for (Int_t r=0; r<yBins; r++)
807 for (Int_t u=0; u<xBins; u++)
808 for (Int_t s=0; s<yBins; s++)
811 covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
812 * (1.0 - hResponse->GetBinContent(u+1, r+1));
814 covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
815 * hResponse->GetBinContent(u+1, s+1);
818 for (Int_t k=0; k<xBins; k++)
819 for (Int_t l=0; l<yBins; l++)
822 printf("In Cov2: %d %d\n", k, l);
823 for (Int_t i=0; i<yBins; i++)
824 for (Int_t j=0; j<yBins; j++)
826 //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
827 // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
829 for (Int_t r=0; r<yBins; r++)
830 for (Int_t u=0; u<xBins; u++)
831 for (Int_t s=0; s<yBins; s++)
833 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
834 || hResponse->GetBinContent(u+1, i+1) == 0)
837 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
838 * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
842 cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
846 printf("Cov2 finished\n");
848 for (Int_t i=1; i<=xBins; i++)
849 for (Int_t j=1; j<=yBins; j++)
850 cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
856 //____________________________________________________________________
857 Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
859 // homogenity term for minuit fitting
860 // pure function of the parameters
861 // prefers constant function (pol0)
863 // Does not take into account efficiency
866 for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
868 Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
869 Double_t left = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
873 Double_t diff = (right - left);
874 chi2 += diff * diff / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
881 //____________________________________________________________________
882 Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
884 // homogenity term for minuit fitting
885 // pure function of the parameters
886 // prefers linear function (pol1)
888 // Does not take into account efficiency
891 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
893 if (params[i-1] == 0)
896 Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
897 Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
898 Double_t left = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
900 Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
901 Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
903 //Double_t diff = (der1 - der2) / middle;
904 //chi2 += diff * diff;
905 chi2 += (der1 - der2) * (der1 - der2) / middle * fgUnfoldedAxis->GetBinWidth(i);
911 //____________________________________________________________________
912 Double_t AliUnfolding::RegularizationLog(TVectorD& params)
914 // homogenity term for minuit fitting
915 // pure function of the parameters
916 // prefers logarithmic function (log)
918 // Does not take into account efficiency
922 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
924 if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
927 Double_t right = log(params[i] / fgUnfoldedAxis->GetBinWidth(i+1));
928 Double_t middle = log(params[i-1] / fgUnfoldedAxis->GetBinWidth(i));
929 Double_t left = log(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1));
931 Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
932 Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
934 //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
936 //if (fgCallCount == 0)
937 // Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
938 chi2 += (der1 - der2) * (der1 - der2);// / error;
944 //____________________________________________________________________
945 Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
947 // homogenity term for minuit fitting
948 // pure function of the parameters
949 // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
950 // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
952 // Does not take into account efficiency
956 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
958 Double_t right = params[i];
959 Double_t middle = params[i-1];
960 Double_t left = params[i-2];
962 Double_t der1 = (right - middle);
963 Double_t der2 = (middle - left);
965 Double_t diff = (der1 - der2);
973 //____________________________________________________________________
974 Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
976 // homogenity term for minuit fitting
977 // pure function of the parameters
978 // calculates entropy, from
979 // The method of reduced cross-entropy (M. Schmelling 1993)
981 // Does not take into account efficiency
983 Double_t paramSum = 0;
985 for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
986 paramSum += params[i];
989 for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
991 Double_t tmp = params[i] / paramSum;
992 //Double_t tmp = params[i];
993 if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
995 chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
1004 //____________________________________________________________________
1005 Double_t AliUnfolding::RegularizationRatio(TVectorD& params)
1007 // homogenity term for minuit fitting
1008 // pure function of the parameters
1010 // Does not take into account efficiency
1014 for (Int_t i=5+fgSkipBinsBegin; i<fgMaxParams; ++i)
1016 if (params[i-1] == 0 || params[i] == 0)
1019 Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1020 Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
1021 Double_t left = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
1022 Double_t left2 = params[i-3] / fgUnfoldedAxis->GetBinWidth(i-2);
1023 Double_t left3 = params[i-4] / fgUnfoldedAxis->GetBinWidth(i-3);
1024 Double_t left4 = params[i-5] / fgUnfoldedAxis->GetBinWidth(i-4);
1026 //Double_t diff = left / middle - middle / right;
1027 //Double_t diff = 2 * left / middle - middle / right - left2 / left;
1028 Double_t diff = 4 * left2 / left - middle / right - left / middle - left3 / left2 - left4 / left3;
1030 chi2 += diff * diff;// / middle;
1036 //____________________________________________________________________
1037 Double_t AliUnfolding::RegularizationPowerLaw(TVectorD& params)
1039 // homogenity term for minuit fitting
1040 // pure function of the parameters
1041 // prefers power law with n = -5
1043 // Does not take into account efficiency
1047 Double_t right = 0.;
1048 Double_t middle = 0.;
1051 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
1053 if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1056 if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i-1) < 1e-8)
1059 middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern);
1062 right = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i),fgPowern)/middle;
1064 left = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle;
1068 Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1069 Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i-2)) / 2);
1071 chi2 += (der1 - der2) * (der1 - der2)/ ( fgUnfoldedAxis->GetBinWidth(i)/2. + fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i-2)/2.)/( fgUnfoldedAxis->GetBinWidth(i)/2. + fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i-2)/2.);// / error;
1072 // printf("i: %d chi2 = %f\n",i,chi2);
1080 //____________________________________________________________________
1081 Double_t AliUnfolding::RegularizationLogLog(TVectorD& params)
1083 // homogenity term for minuit fitting
1084 // pure function of the parameters
1085 // prefers a powerlaw (linear on a log-log scale)
1087 // The calculation takes into account the efficiencies
1091 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
1093 if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
1095 if ((*fgEfficiency)(i-1) == 0 || (*fgEfficiency)(i) == 0 || (*fgEfficiency)(i-2) == 0)
1099 Double_t right = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i));
1100 Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i-1));
1101 Double_t left = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-2));
1103 Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) );
1104 Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) );
1106 double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.;
1107 Double_t dder = (der1-der2) / tmp;
1109 chi2 += dder * dder;
1117 //____________________________________________________________________
1118 void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
1121 // fit function for minuit
1122 // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
1125 // TODO use static members for the variables here to speed up processing (no construction/deconstruction)
1128 TVectorD paramsVector(fgMaxParams);
1129 for (Int_t i=0; i<fgMaxParams; ++i)
1130 paramsVector[i] = params[i] * params[i];
1132 // calculate penalty factor
1133 Double_t penaltyVal = 0;
1135 switch (fgRegularizationType)
1138 case kPol0: penaltyVal = RegularizationPol0(paramsVector); break;
1139 case kPol1: penaltyVal = RegularizationPol1(paramsVector); break;
1140 case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
1141 case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
1142 case kLog: penaltyVal = RegularizationLog(paramsVector); break;
1143 case kRatio: penaltyVal = RegularizationRatio(paramsVector); break;
1144 case kPowerLaw: penaltyVal = RegularizationPowerLaw(paramsVector); break;
1145 case kLogLog: penaltyVal = RegularizationLogLog(paramsVector); break;
1148 penaltyVal *= fgRegularizationWeight; //beta*PU
1151 TVectorD measGuessVector(fgMaxInput);
1152 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1155 measGuessVector -= (*fgCurrentESDVector);
1158 // new error calcuation using error on the guess
1161 TVectorD errorGuessVector(fgMaxInput);
1162 //errorGuessVector = (*fgCorrelationMatrixSquared) * paramsVector;
1163 errorGuessVector = (*fgCorrelationMatrix) * paramsVector;
1165 Double_t chi2FromFit = 0;
1166 for (Int_t i=0; i<fgMaxInput; ++i)
1169 if (errorGuessVector(i) > 0)
1170 error = errorGuessVector(i);
1171 chi2FromFit += measGuessVector(i) * measGuessVector(i) / error;
1175 // old error calcuation using the error on the measured
1176 TVectorD copy(measGuessVector);
1179 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1180 // normal way is like this:
1181 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1182 // optimized way like this:
1183 for (Int_t i=0; i<fgMaxInput; ++i)
1184 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1187 if (fgSkipBin0InChi2)
1188 measGuessVector[0] = 0;
1190 // (Ad - m) W (Ad - m)
1191 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1192 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see UnfoldWithMinuit)
1193 Double_t chi2FromFit = measGuessVector * copy * 1e6;
1196 Double_t notFoundEventsConstraint = 0;
1197 Double_t currentNotFoundEvents = 0;
1198 Double_t errorNotFoundEvents = 0;
1200 if (fgNotFoundEvents > 0)
1202 for (Int_t i=0; i<fgMaxParams; ++i)
1204 Double_t eff = (1.0 / (*fgEfficiency)(i) - 1);
1206 eff = (1.0 / params[fgMaxParams] - 1);
1207 currentNotFoundEvents += eff * paramsVector(i);
1208 errorNotFoundEvents += eff * eff * paramsVector(i); // error due to guess (paramsVector)
1210 errorNotFoundEvents += (eff * 0.03) * (eff * 0.03) * paramsVector(i) * paramsVector(i); // error on eff
1211 // if ((fgCallCount % 10000) == 0)
1212 //Printf("%d %f %f %f", i, (*fgEfficiency)(i), paramsVector(i), currentNotFoundEvents);
1214 errorNotFoundEvents += fgNotFoundEvents;
1215 // TODO add error on background, background estimate
1217 notFoundEventsConstraint = (currentNotFoundEvents - fgNotFoundEvents) * (currentNotFoundEvents - fgNotFoundEvents) / errorNotFoundEvents;
1222 Float_t currentMult = 0;
1223 // try to match dn/deta
1224 for (Int_t i=0; i<fgMaxParams; ++i)
1226 avg += paramsVector(i) * currentMult;
1227 sum += paramsVector(i);
1228 currentMult += fgUnfoldedAxis->GetBinWidth(i);
1231 Float_t chi2Avg = 0; //(avg - 3.73) * (avg - 3.73) * 100;
1233 chi2 = chi2FromFit + penaltyVal + notFoundEventsConstraint + chi2Avg;
1235 if ((fgCallCount++ % 1000) == 0)
1238 Printf("AliUnfolding::Chi2Function: Iteration %d (ev %d %d +- %f) (%f) (%f): %f %f %f %f --> %f", fgCallCount-1, (Int_t) fgNotFoundEvents, (Int_t) currentNotFoundEvents, TMath::Sqrt(errorNotFoundEvents), params[fgMaxParams-1], avg, chi2FromFit, penaltyVal, notFoundEventsConstraint, chi2Avg, chi2);
1240 //for (Int_t i=0; i<fgMaxInput; ++i)
1241 // Printf("%d: %f", i, measGuessVector(i) * copy(i) * 1e6);
1244 fChi2FromFit = chi2FromFit;
1245 fPenaltyVal = penaltyVal;
1248 //____________________________________________________________________
1249 void AliUnfolding::MakePads() {
1250 TPad *presult = new TPad("presult","result",0,0.4,1,1);
1251 presult->SetNumber(1);
1254 TPad *pres = new TPad("pres","residuals",0,0.2,1,0.4);
1257 TPad *ppen = new TPad("ppen","penalty",0,0.0,1,0.2);
1262 //____________________________________________________________________
1263 void AliUnfolding::DrawResults(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TCanvas *canv, Int_t reuseHists,TH1 *unfolded)
1265 // Draw histograms of
1266 // - Result folded with response
1267 // - Penalty factors
1268 // - Chisquare contributions
1269 // (Useful for debugging/sanity checks and the interactive unfolder)
1271 // If a canvas pointer is given (canv != 0), it will be used for all
1272 // plots; 3 pads are made if needed.
1275 Int_t blankCanvas = 0;
1276 TVirtualPad *presult = 0;
1277 TVirtualPad *pres = 0;
1278 TVirtualPad *ppen = 0;
1282 presult = canv->GetPad(1);
1283 pres = canv->GetPad(2);
1284 ppen = canv->GetPad(3);
1285 if (presult == 0 || pres == 0 || ppen == 0) {
1289 presult = canv->GetPad(1);
1290 pres = canv->GetPad(2);
1291 ppen = canv->GetPad(3);
1295 presult = new TCanvas;
1301 if (fgMaxInput == -1)
1303 Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
1304 fgMaxInput = measured->GetNbinsX();
1306 if (fgMaxParams == -1)
1308 Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
1309 fgMaxParams = initialConditions->GetNbinsX();
1312 if (fgOverflowBinLimit > 0)
1313 CreateOverflowBin(correlation, measured);
1315 // Uses Minuit methods
1317 SetStaticVariables(correlation, measured, efficiency);
1319 Double_t* params = new Double_t[fgMaxParams+1];
1320 for (Int_t i=0; i<fgMaxParams; ++i)
1322 params[i] = initialConditions->GetBinContent(i+1) * efficiency->GetBinContent(i+1);
1324 Bool_t fix = kFALSE;
1328 params[i] = -params[i];
1330 params[i]=TMath::Sqrt(params[i]);
1332 //cout << "params[" << i << "] " << params[i] << endl;
1339 //fgPrintChi2Details = kTRUE;
1340 fgCallCount = 0; // To make sure that Chi2Function prints the components
1341 Chi2Function(dummy, 0, chi2, params, 0);
1344 TH1 *meas2 = (TH1*)measured->Clone("meas2");
1345 meas2->SetTitle("meas2");
1346 meas2->SetName("meas2");
1347 meas2->SetLineColor(1);
1348 meas2->SetLineWidth(2);
1349 meas2->SetMarkerColor(meas2->GetLineColor());
1350 meas2->SetMarkerStyle(20);
1351 Float_t scale = unfolded->GetXaxis()->GetBinWidth(1)/meas2->GetXaxis()->GetBinWidth(1);
1352 meas2->Scale(scale);
1356 meas2->DrawCopy("same");
1357 //meas2->GetXaxis()->SetLimits(0,fgMaxInput);
1358 meas2->SetBit(kCannotPick);
1359 DrawGuess(params, presult, pres, ppen, reuseHists,unfolded);
1362 //____________________________________________________________________
1363 void AliUnfolding::RedrawInteractive() {
1365 // Helper function for interactive unfolding
1367 DrawResults(fghCorrelation,fghEfficiency,fghMeasured,fghUnfolded,fgCanvas,1,fghUnfolded);
1369 //____________________________________________________________________
1370 void AliUnfolding::InteractiveUnfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions)
1373 // Function to do interactive unfolding
1374 // A canvas is drawn with the unfolding result
1375 // Change the histogram with your mouse and all histograms
1376 // will be updated automatically
1378 fgCanvas = new TCanvas("UnfoldingCanvas","Interactive unfolding",500,800);
1384 delete fghUnfolded; fghUnfolded = 0;
1387 gROOT->SetEditHistograms(kTRUE);
1389 fghUnfolded = new TH1F("AliUnfoldingfghUnfolded","Unfolded result (Interactive unfolder",efficiency->GetNbinsX(),efficiency->GetXaxis()->GetXmin(),efficiency->GetXaxis()->GetXmax());
1391 fghCorrelation = correlation;
1392 fghEfficiency = efficiency;
1393 fghMeasured = measured;
1395 if(fgMinuitMaxIterations>0)
1396 UnfoldWithMinuit(correlation,efficiency,measured,initialConditions,fghUnfolded,kFALSE);
1398 fghUnfolded = initialConditions;
1400 fghUnfolded->SetLineColor(2);
1401 fghUnfolded->SetMarkerColor(2);
1402 fghUnfolded->SetLineWidth(2);
1406 fghUnfolded->Draw();
1407 DrawResults(correlation,efficiency,measured,fghUnfolded,fgCanvas,kFALSE,fghUnfolded);
1409 TExec *execRedraw = new TExec("redraw","AliUnfolding::RedrawInteractive()");
1410 fghUnfolded->GetListOfFunctions()->Add(execRedraw);
1412 //____________________________________________________________________
1413 void AliUnfolding::DrawGuess(Double_t *params, TVirtualPad *pfolded, TVirtualPad *pres, TVirtualPad *ppen, Int_t reuseHists,TH1* unfolded)
1416 // draws residuals of solution suggested by params and effect of regularization
1420 pfolded = new TCanvas;
1427 TVectorD paramsVector(fgMaxParams);
1428 for (Int_t i=0; i<fgMaxParams; ++i)
1429 paramsVector[i] = params[i] * params[i];
1432 TVectorD measGuessVector(fgMaxInput);
1433 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1437 folded = dynamic_cast<TH1*>(gROOT->FindObject("__hfolded"));
1438 if (!reuseHists || folded == 0) {
1439 if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins
1440 folded = new TH1F("__hfolded","Folded histo from AliUnfolding",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray());
1442 folded = new TH1F("__hfolded","Folded histo from AliUnfolding",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(),fgMeasuredAxis->GetXmax());
1445 folded->SetBit(kCannotPick);
1446 folded->SetLineColor(4);
1447 folded->SetLineWidth(2);
1449 for (Int_t ibin =0; ibin < fgMaxInput; ibin++)
1450 folded->SetBinContent(ibin+1, measGuessVector[ibin]);
1452 Float_t scale = unfolded->GetXaxis()->GetBinWidth(1)/folded->GetXaxis()->GetBinWidth(1);
1453 folded->Scale(scale);
1457 folded->Draw("same");
1460 measGuessVector -= (*fgCurrentESDVector);
1462 TVectorD copy(measGuessVector);
1465 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1466 // normal way is like this:
1467 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1468 // optimized way like this:
1469 for (Int_t i=0; i<fgMaxInput; ++i)
1470 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1472 // (Ad - m) W (Ad - m)
1473 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1474 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
1475 //Double_t chi2FromFit = measGuessVector * copy * 1e6;
1478 // Double_t pTarray[fgMaxParams+1];
1479 // for(int i=0; i<fgMaxParams; i++) {
1480 // pTarray[i] = fgUnfoldedAxis->GetBinCenter(i)-0.5*fgUnfoldedAxis->GetBinWidth(i);
1482 // pTarray[fgMaxParams] = fgUnfoldedAxis->GetBinCenter(fgMaxParams-1)+0.5*fgUnfoldedAxis->GetBinWidth(fgMaxParams-1);
1483 // TH1* residuals = new TH1F("residuals", "residuals", fgMaxParams,pTarray);
1484 // TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
1485 // for (Int_t i=0; i<fgMaxInput; ++i)
1486 // residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);7
1487 TH1* residuals = GetResidualsPlot(params);
1488 residuals->GetXaxis()->SetTitleSize(0.06);
1489 residuals->GetXaxis()->SetTitleOffset(0.7);
1490 residuals->GetXaxis()->SetLabelSize(0.07);
1491 residuals->GetYaxis()->SetTitleSize(0.08);
1492 residuals->GetYaxis()->SetTitleOffset(0.5);
1493 residuals->GetYaxis()->SetLabelSize(0.06);
1494 pres->cd(); residuals->DrawCopy(); gPad->SetLogy();
1498 TH1* penalty = GetPenaltyPlot(params);
1499 penalty->GetXaxis()->SetTitleSize(0.06);
1500 penalty->GetXaxis()->SetTitleOffset(0.7);
1501 penalty->GetXaxis()->SetLabelSize(0.07);
1502 penalty->GetYaxis()->SetTitleSize(0.08);
1503 penalty->GetYaxis()->SetTitleOffset(0.5);
1504 penalty->GetYaxis()->SetLabelSize(0.06);
1505 ppen->cd(); penalty->DrawCopy(); gPad->SetLogy();
1508 //____________________________________________________________________
1509 TH1* AliUnfolding::GetResidualsPlot(TH1* corrected)
1512 // MvL: THIS MUST BE INCORRECT.
1513 // Need heff to calculate params from TH1 'corrected'
1517 // fill residuals histogram of solution suggested by params and effect of regularization
1520 Double_t* params = new Double_t[fgMaxParams];
1521 for (Int_t i=0; i<fgMaxParams; i++)
1524 for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
1525 params[i] = TMath::Sqrt(TMath::Abs(corrected->GetBinContent(i+1)*(*fgEfficiency)(i)));
1527 TH1 * plot = GetResidualsPlot(params);
1532 //____________________________________________________________________
1533 TH1* AliUnfolding::GetResidualsPlot(Double_t* params)
1536 // fill residuals histogram of solution suggested by params and effect of regularization
1540 TVectorD paramsVector(fgMaxParams);
1541 for (Int_t i=0; i<fgMaxParams; ++i)
1542 paramsVector[i] = params[i] * params[i];
1545 TVectorD measGuessVector(fgMaxInput);
1546 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1549 measGuessVector -= (*fgCurrentESDVector);
1551 TVectorD copy(measGuessVector);
1554 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1555 // normal way is like this:
1556 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1557 // optimized way like this:
1558 for (Int_t i=0; i<fgMaxInput; ++i)
1559 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1561 // (Ad - m) W (Ad - m)
1562 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1563 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
1564 //Double_t chi2FromFit = measGuessVector * copy * 1e6;
1568 if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins
1569 residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray());
1571 residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(), fgMeasuredAxis->GetXmax());
1572 // TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
1574 Double_t sumResiduals = 0.;
1575 for (Int_t i=0; i<fgMaxInput; ++i) {
1576 residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
1577 sumResiduals += measGuessVector(i) * copy(i) * 1e6;
1579 fAvgResidual = sumResiduals/(double)fgMaxInput;
1584 //____________________________________________________________________
1585 TH1* AliUnfolding::GetPenaltyPlot(TH1* corrected)
1587 // draws the penalty factors as function of multiplicity of the current selected regularization
1589 Double_t* params = new Double_t[fgMaxParams];
1590 for (Int_t i=0; i<fgMaxParams; i++)
1593 for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
1594 params[i] = (*fgEfficiency)(i)*corrected->GetBinContent(i+1);
1596 TH1* penalty = GetPenaltyPlot(params);
1603 //____________________________________________________________________
1604 TH1* AliUnfolding::GetPenaltyPlot(Double_t* params)
1606 // draws the penalty factors as function of multiplicity of the current selected regularization
1608 //TH1* penalty = new TH1F("penalty", ";unfolded multiplicity;penalty factor", fgMaxParams, -0.5, fgMaxParams - 0.5);
1609 // TH1* penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgMaxParams, fgUnfoldedAxis->GetBinCenter(0)-0.5*fgUnfoldedAxis->GetBinWidth(0),fgUnfoldedAxis->GetBinCenter(fgMaxParams)+0.5*fgUnfoldedAxis->GetBinWidth(fgMaxParams) );
1612 if (fgUnfoldedAxis->GetXbins()->GetArray())
1613 penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXbins()->GetArray());
1615 penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXmin(),fgUnfoldedAxis->GetXmax());
1617 for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
1620 if (fgRegularizationType == kPol0)
1622 Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1623 Double_t left = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
1627 Double_t diffTmp = (right - left);
1628 diff = diffTmp * diffTmp / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2) / 100;
1631 if (fgRegularizationType == kPol1 && i > 1+fgSkipBinsBegin)
1633 if (params[i-1] == 0)
1636 Double_t right = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1637 Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
1638 Double_t left = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
1640 Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1641 Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
1643 diff = (der1 - der2) * (der1 - der2) / middle;
1646 if (fgRegularizationType == kLog && i > 1+fgSkipBinsBegin)
1648 if (params[i-1] == 0)
1651 Double_t right = log(params[i]);
1652 Double_t middle = log(params[i-1]);
1653 Double_t left = log(params[i-2]);
1655 Double_t der1 = (right - middle);
1656 Double_t der2 = (middle - left);
1658 //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
1659 //Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
1661 diff = (der1 - der2) * (der1 - der2);// / error;
1663 if (fgRegularizationType == kCurvature && i > 1+fgSkipBinsBegin)
1665 Double_t right = params[i]; // params are sqrt
1666 Double_t middle = params[i-1];
1667 Double_t left = params[i-2];
1669 Double_t der1 = (right - middle)/0.5/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i));
1670 Double_t der2 = (middle - left)/0.5/(fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i+1));
1672 diff = (der1 - der2)/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1))*3.0;
1673 diff = 1e4*diff*diff;
1675 if (fgRegularizationType == kPowerLaw && i > 1+fgSkipBinsBegin)
1678 if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1681 if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8)
1684 double middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern);
1687 double right = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i+1),fgPowern)/middle;
1689 double left = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle;
1693 Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1694 Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
1696 diff = (der1 - der2) * (der1 - der2);// / error;
1700 if (fgRegularizationType == kLogLog && i > 1+fgSkipBinsBegin)
1703 if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1706 Double_t right = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i+1));
1707 Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i));
1708 Double_t left = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-1));
1710 Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) );
1711 Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) );
1713 double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.;
1714 Double_t dder = (der1-der2) / tmp;
1719 penalty->SetBinContent(i, diff*fgRegularizationWeight);
1721 //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
1727 //____________________________________________________________________
1728 void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
1731 // fit function for minuit
1732 // uses the TF1 stored in fgFitFunction
1735 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1736 fgFitFunction->SetParameter(i, params[i]);
1738 Double_t* params2 = new Double_t[fgMaxParams];
1740 for (Int_t i=0; i<fgMaxParams; ++i)
1741 params2[i] = fgFitFunction->Eval(i);
1743 Chi2Function(unused1, unused2, chi2, params2, unused3);
1751 //____________________________________________________________________
1752 Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult)
1755 // correct spectrum using minuit chi2 method applying a functional fit
1760 Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting.");
1764 SetChi2Regularization(kNone, 0);
1766 SetStaticVariables(correlation, measured, efficiency);
1768 // Initialize TMinuit via generic fitter interface
1769 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
1771 minuit->SetFCN(TF1Function);
1772 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1774 Double_t lower, upper;
1775 fgFitFunction->GetParLimits(i, lower, upper);
1776 minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper);
1779 Double_t arglist[100];
1781 minuit->ExecuteCommand("SET PRINT", arglist, 1);
1782 minuit->ExecuteCommand("SCAN", arglist, 0);
1783 minuit->ExecuteCommand("MIGRAD", arglist, 0);
1784 //minuit->ExecuteCommand("MINOS", arglist, 0);
1786 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1787 fgFitFunction->SetParameter(i, minuit->GetParameter(i));
1789 for (Int_t i=0; i<fgMaxParams; ++i)
1791 Double_t value = fgFitFunction->Eval(i);
1793 Printf("%d : %f", i, value);
1797 if (efficiency->GetBinContent(i+1) > 0)
1799 value /= efficiency->GetBinContent(i+1);
1804 aResult->SetBinContent(i+1, value);
1805 aResult->SetBinError(i+1, 0);
1811 //____________________________________________________________________
1812 void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
1814 // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins
1815 // The same limit is applied to the correlation
1818 for (Int_t i=1; i<=measured->GetNbinsX(); ++i)
1820 if (measured->GetBinContent(i) <= fgOverflowBinLimit)
1829 Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit);
1833 Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin);
1835 TCanvas* canvas = 0;
1839 canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
1840 canvas->Divide(2, 2);
1843 measured->SetStats(kFALSE);
1844 measured->DrawCopy();
1848 correlation->SetStats(kFALSE);
1849 correlation->DrawCopy("COLZ");
1852 measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX()));
1853 for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i)
1855 measured->SetBinContent(i, 0);
1856 measured->SetBinError(i, 0);
1858 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1859 measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin)));
1861 Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin));
1863 for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1865 correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY()));
1866 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1867 correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin)));
1869 for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j)
1871 correlation->SetBinContent(i, j, 0);
1872 correlation->SetBinError(i, j, 0);
1876 Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!");
1881 measured->DrawCopy();
1885 correlation->DrawCopy("COLZ");
1889 Int_t AliUnfolding::UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result)
1891 // unfolds and assigns bias as errors with Eq. 19 of Cowan, "a survey of unfolding methods for particle physics"
1892 // b_i = sum_j dmu_i/dn_j (nu_j - n_j) with nu_j as folded guess, n_j as data
1893 // dmu_i/dn_j is found numerically by changing the bin content and re-unfolding
1895 // return code: 0 (success) -1 (error: from Unfold(...) )
1897 if (Unfold(correlation, efficiency, measured, initialConditions, result) != 0)
1900 TMatrixD matrix(fgMaxInput, fgMaxParams);
1903 for (Int_t loop=0; loop<4; loop++)
1904 newResult[loop] = (TH1*) result->Clone(Form("newresult_%d", loop));
1906 // change bin-by-bin and built matrix of effects
1907 for (Int_t m=0; m<fgMaxInput; m++)
1909 if (measured->GetBinContent(m+1) < 1)
1912 for (Int_t loop=0; loop<4; loop++)
1914 //Printf("%d %d", i, loop);
1918 case 0: factor = 0.5; break;
1919 case 1: factor = -0.5; break;
1920 case 2: factor = 1; break;
1921 case 3: factor = -1; break;
1925 TH1* measuredClone = (TH1*) measured->Clone("measuredClone");
1927 measuredClone->SetBinContent(m+1, measured->GetBinContent(m+1) + factor * TMath::Sqrt(measured->GetBinContent(m+1)));
1928 //new TCanvas; measuredClone->Draw("COLZ");
1930 newResult[loop]->Reset();
1931 Unfold(correlation, efficiency, measuredClone, measuredClone, newResult[loop]);
1932 // WARNING if we leave here, then nothing is calculated
1935 delete measuredClone;
1938 for (Int_t t=0; t<fgMaxParams; t++)
1941 //matrix(i, j) = (result->GetBinContent(j+1) - newResult->GetBinContent(j+1)) / TMath::Sqrt(changedHist->GetBinContent(1, i+1));
1943 // four values from the derivate (procedure taken from ROOT -- suggestion by Ruben)
1944 // = 1/2D * [ 8 (f(D/2) - f(-D/2)) - (f(D)-f(-D)) ]/3
1948 measured->SetBinContent(1, m+1, 100);
1949 newResult[0]->SetBinContent(t+1, 5 * 0.5 + 10);
1950 newResult[1]->SetBinContent(t+1, 5 * -0.5 + 10);
1951 newResult[2]->SetBinContent(t+1, 5 * 1 + 10);
1952 newResult[3]->SetBinContent(t+1, 5 * -1 + 10);
1955 matrix(m, t) = 0.5 / TMath::Sqrt(measured->GetBinContent(m+1)) *
1956 ( 8. * (newResult[0]->GetBinContent(t+1) - newResult[1]->GetBinContent(t+1)) -
1957 (newResult[2]->GetBinContent(t+1) - newResult[3]->GetBinContent(t+1))
1962 for (Int_t loop=0; loop<4; loop++)
1963 delete newResult[loop];
1965 // ...calculate folded guess
1966 TH1* convoluted = (TH1*) measured->Clone("convoluted");
1967 convoluted->Reset();
1968 for (Int_t m=0; m<fgMaxInput; m++)
1971 for (Int_t t = 0; t<fgMaxParams; t++)
1973 Float_t tmp = correlation->GetBinContent(t+1, m+1) * result->GetBinContent(t+1);
1975 tmp *= efficiency->GetBinContent(t+1);
1978 convoluted->SetBinContent(m+1, value);
1982 convoluted->Scale(measured->Integral() / convoluted->Integral());
1984 //new TCanvas; convoluted->Draw(); measured->Draw("SAME"); measured->SetLineColor(2);
1988 convoluted->Add(measured, -1);
1991 for (Int_t t = 0; t<fgMaxParams; t++)
1994 for (Int_t m=0; m<fgMaxInput; m++)
1995 error += matrix(m, t) * convoluted->GetBinContent(m+1);
1996 result->SetBinError(t+1, error);
1999 //new TCanvas; result->Draw(); gPad->SetLogy();