Added check of daughters clustermap in ITS
[u/mrichter/AliRoot.git] / PWG0 / AliUnfolding.cxx
CommitLineData
19442b86 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16/* $Id: AliUnfolding.cxx 31168 2009-02-23 15:18:45Z jgrosseo $ */
17
18// This class allows 1-dimensional unfolding.
19// Methods that are implemented are chi2 minimization and bayesian unfolding.
20//
21// Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
22
23#include "AliUnfolding.h"
24#include <TH1F.h>
25#include <TH2F.h>
26#include <TVirtualFitter.h>
27#include <TMath.h>
28#include <TCanvas.h>
29#include <TF1.h>
30
31TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
32TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
33TVectorD* AliUnfolding::fgCurrentESDVector = 0;
34TVectorD* AliUnfolding::fgEntropyAPriori = 0;
35
36TF1* AliUnfolding::fgFitFunction = 0;
37
38AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid;
39Int_t AliUnfolding::fgMaxInput = -1; // bins in measured histogram
40Int_t AliUnfolding::fgMaxParams = -1; // bins in unfolded histogram = number of fit params
41Float_t AliUnfolding::fgOverflowBinLimit = -1;
42
43AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1;
44Float_t AliUnfolding::fgRegularizationWeight = 10000;
45Int_t AliUnfolding::fgSkipBinsBegin = 0;
46Float_t AliUnfolding::fgMinuitStepSize = 0.1; // (usually not needed to be changed) step size in minimization
47Bool_t AliUnfolding::fgNormalizeInput = kFALSE; // normalize input spectrum
48
49Float_t AliUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
651e2088 50Int_t AliUnfolding::fgBayesianIterations = 10; // number of iterations in Bayesian method
19442b86 51
52Bool_t AliUnfolding::fgDebug = kFALSE;
53
54ClassImp(AliUnfolding)
55
56//____________________________________________________________________
57void AliUnfolding::SetUnfoldingMethod(MethodType methodType)
58{
59 // set unfolding method
60 fgMethodType = methodType;
61
62 const char* name = 0;
63 switch (methodType)
64 {
65 case kInvalid: name = "INVALID"; break;
66 case kChi2Minimization: name = "Chi2 Minimization"; break;
67 case kBayesian: name = "Bayesian unfolding"; break;
68 case kFunction: name = "Functional fit"; break;
69 }
70 Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name);
71}
72
73//____________________________________________________________________
74void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit)
75{
76 // enable the creation of a overflow bin that includes all statistics below the given limit
77
78 fgOverflowBinLimit = overflowBinLimit;
79
80 Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit);
81}
82
83//____________________________________________________________________
84void AliUnfolding::SetSkipBinsBegin(Int_t nBins)
85{
86 // set number of skipped bins in regularization
87
88 fgSkipBinsBegin = nBins;
89
90 Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin);
91}
92
93//____________________________________________________________________
94void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded)
95{
96 // set number of bins in the input (measured) distribution and in the unfolded distribution
97 fgMaxInput = nMeasured;
98 fgMaxParams = nUnfolded;
99
100 Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
101}
102
103//____________________________________________________________________
104void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight)
105{
106 //
107 // sets the parameters for chi2 minimization
108 //
109
110 fgRegularizationType = type;
111 fgRegularizationWeight = weight;
112
113 Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight);
114}
115
116//____________________________________________________________________
117void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
118{
119 //
120 // sets the parameters for Bayesian unfolding
121 //
122
123 fgBayesianSmoothing = smoothing;
124 fgBayesianIterations = nIterations;
125
126 Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
127}
128
129//____________________________________________________________________
130void AliUnfolding::SetFunction(TF1* function)
131{
132 // set function for unfolding with a fit function
133
134 fgFitFunction = function;
135
136 Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar());
137}
138
139//____________________________________________________________________
140Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
141{
142 // unfolds with unfolding method fgMethodType
143 //
144 // parameters:
145 // correlation: response matrix as measured vs. generated
146 // 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.
147 // measured: the measured spectrum
148 // initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions.
149 // result: target for the unfolded result
150 // check: depends on the unfolding method, see comments in specific functions
151
152 if (fgMaxInput == -1)
153 {
154 Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
155 fgMaxInput = measured->GetNbinsX();
156 }
157 if (fgMaxParams == -1)
158 {
159 Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
160 fgMaxParams = measured->GetNbinsX();
161 }
162
163 if (fgOverflowBinLimit > 0)
164 CreateOverflowBin(correlation, measured);
165
166 switch (fgMethodType)
167 {
168 case kInvalid:
169 {
170 Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting...");
171 return -1;
172 }
173 case kChi2Minimization:
174 return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check);
175 case kBayesian:
176 return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result);
177 case kFunction:
178 return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result);
179 }
180 return -1;
181}
182
183//____________________________________________________________________
184void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured)
185{
186 // fill static variables needed for minuit fit
187
188 if (!fgCorrelationMatrix)
189 fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
190 if (!fgCorrelationCovarianceMatrix)
191 fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
192 if (!fgCurrentESDVector)
193 fgCurrentESDVector = new TVectorD(fgMaxInput);
194 if (!fgEntropyAPriori)
195 fgEntropyAPriori = new TVectorD(fgMaxParams);
196
197 fgCorrelationMatrix->Zero();
198 fgCorrelationCovarianceMatrix->Zero();
199 fgCurrentESDVector->Zero();
200 fgEntropyAPriori->Zero();
201
202 // normalize correction for given nPart
203 for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
204 {
205 Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
206 if (sum <= 0)
207 continue;
208 Float_t maxValue = 0;
209 Int_t maxBin = -1;
210 for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
211 {
212 // find most probably value
213 if (maxValue < correlation->GetBinContent(i, j))
214 {
215 maxValue = correlation->GetBinContent(i, j);
216 maxBin = j;
217 }
218
219 // npart sum to 1
220 correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
221 correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
222
223 if (i <= fgMaxParams && j <= fgMaxInput)
224 (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
225 }
226
227 //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
228 }
229
230 //normalize measured
231 if (fgNormalizeInput)
232 measured->Scale(1.0 / measured->Integral());
233
234 for (Int_t i=0; i<fgMaxInput; ++i)
235 {
236 (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
237 if (measured->GetBinError(i+1) > 0)
238 (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
239
240 if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
241 (*fgCorrelationCovarianceMatrix)(i, i) = 0;
242 //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
243 }
244}
245
246//____________________________________________________________________
247Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
248{
249 //
250 // implementation of unfolding (internal function)
251 //
252 // unfolds <measured> using response from <correlation> and effiency <efficiency>
253 // output is in <result>
254 // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
255 // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
256 //
257 // returns minuit status (0 = success), (-1 when check was set)
258 //
259
260 SetStaticVariables(correlation, measured);
261
262 // Initialize TMinuit via generic fitter interface
263 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
264 Double_t arglist[100];
265
266 // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
267 arglist[0] = 0;
268 minuit->ExecuteCommand("SET PRINT", arglist, 1);
269
270 // however, enable warnings
271 //minuit->ExecuteCommand("SET WAR", arglist, 0);
272
273 // set minimization function
274 minuit->SetFCN(Chi2Function);
275
276 for (Int_t i=0; i<fgMaxParams; i++)
277 (*fgEntropyAPriori)[i] = 1;
278
279 if (!initialConditions) {
280 initialConditions = measured;
281 } else {
282 Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions...");
283 //new TCanvas; initialConditions->DrawCopy();
284 if (fgNormalizeInput)
285 initialConditions->Scale(1.0 / initialConditions->Integral());
286 }
287
288 // set initial conditions as a-priori distribution for MRX regularization
289 for (Int_t i=0; i<fgMaxParams; i++)
290 if (initialConditions->GetBinContent(i+1) > 0)
291 (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
292
293 Double_t* results = new Double_t[fgMaxParams+1];
294 for (Int_t i=0; i<fgMaxParams; ++i)
295 {
296 results[i] = initialConditions->GetBinContent(i+1);
297
298 // minuit sees squared values to prevent it from going negative...
299 results[i] = TMath::Sqrt(results[i]);
300
301 minuit->SetParameter(i, Form("param%d", i), results[i], fgMinuitStepSize, 0, 0);
302 }
303
304 Int_t dummy = 0;
305 Double_t chi2 = 0;
306 Chi2Function(dummy, 0, chi2, results, 0);
307 printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2);
308
309 if (check)
310 {
311 DrawGuess(results);
312 delete[] results;
313 return -1;
314 }
315
316 // first param is number of iterations, second is precision....
317 arglist[0] = 1e6;
318 //arglist[1] = 1e-5;
319 //minuit->ExecuteCommand("SCAN", arglist, 0);
320 Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
321 Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status);
322 //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
323 //minuit->ExecuteCommand("MINOS", arglist, 0);
324
325 for (Int_t i=0; i<fgMaxParams; ++i)
326 {
327 results[i] = minuit->GetParameter(i);
328 Double_t value = results[i] * results[i];
329 // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
330 Double_t error = minuit->GetParError(i) * results[i];
331
332 if (efficiency)
333 {
334 if (efficiency->GetBinContent(i+1) > 0)
335 {
336 value /= efficiency->GetBinContent(i+1);
337 error /= efficiency->GetBinContent(i+1);
338 }
339 else
340 {
341 value = 0;
342 error = 0;
343 }
344 }
345
346 result->SetBinContent(i+1, value);
347 result->SetBinError(i+1, error);
348 }
349
350 if (fgDebug)
351 DrawGuess(results);
352
353 delete[] results;
354
355 return status;
356}
357
358//____________________________________________________________________
359Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult)
360{
361 //
362 // unfolds a spectrum using the Bayesian method
363 //
364
365 if (measured->Integral() <= 0)
366 {
367 Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
368 return -1;
369 }
370
371 const Int_t kStartBin = 0;
372
373 Int_t kMaxM = fgMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
374 Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
375
376 // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
377 const Double_t kConvergenceLimit = kMaxT * 1e-6;
378
379 // store information in arrays, to increase processing speed (~ factor 5)
380 Double_t* measuredCopy = new Double_t[kMaxM];
381 Double_t* measuredError = new Double_t[kMaxM];
382 Double_t* prior = new Double_t[kMaxT];
383 Double_t* result = new Double_t[kMaxT];
384 Double_t* efficiency = new Double_t[kMaxT];
385
386 Double_t** response = new Double_t*[kMaxT];
387 Double_t** inverseResponse = new Double_t*[kMaxT];
388 for (Int_t i=0; i<kMaxT; i++)
389 {
390 response[i] = new Double_t[kMaxM];
391 inverseResponse[i] = new Double_t[kMaxM];
392 }
393
394 // for normalization
395 Float_t measuredIntegral = measured->Integral();
396 for (Int_t m=0; m<kMaxM; m++)
397 {
398 measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
399 measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
400
401 for (Int_t t=0; t<kMaxT; t++)
402 {
403 response[t][m] = correlation->GetBinContent(t+1, m+1);
404 inverseResponse[t][m] = 0;
405 }
406 }
407
408 for (Int_t t=0; t<kMaxT; t++)
409 {
651e2088 410 if (aEfficiency)
19442b86 411 {
412 efficiency[t] = aEfficiency->GetBinContent(t+1);
413 }
414 else
415 efficiency[t] = 1;
416
417 prior[t] = measuredCopy[t];
418 result[t] = 0;
419 }
420
421 // pick prior distribution
422 if (initialConditions)
423 {
424 printf("Using different starting conditions...\n");
425 // for normalization
426 Float_t inputDistIntegral = initialConditions->Integral();
427 for (Int_t i=0; i<kMaxT; i++)
428 prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
429 }
430
431 //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
432
433 // unfold...
434 for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations< 0; i++)
435 {
436 if (fgDebug)
437 Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
438
439 // calculate IR from Beyes theorem
440 // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
441
442 Double_t chi2Measured = 0;
443 for (Int_t m=0; m<kMaxM; m++)
444 {
445 Float_t norm = 0;
446 for (Int_t t = kStartBin; t<kMaxT; t++)
447 norm += response[t][m] * prior[t];
448
449 // calc. chi2: (measured - response * prior) / error
450 if (measuredError[m] > 0)
451 {
452 Double_t value = (measuredCopy[m] - norm) / measuredError[m];
453 chi2Measured += value * value;
454 }
455
456 if (norm > 0)
457 {
458 for (Int_t t = kStartBin; t<kMaxT; t++)
459 inverseResponse[t][m] = response[t][m] * prior[t] / norm;
460 }
461 else
462 {
463 for (Int_t t = kStartBin; t<kMaxT; t++)
464 inverseResponse[t][m] = 0;
465 }
466 }
467 //Printf("chi2Measured of the last prior is %e", chi2Measured);
468
469 for (Int_t t = kStartBin; t<kMaxT; t++)
470 {
471 Float_t value = 0;
472 for (Int_t m=0; m<kMaxM; m++)
473 value += inverseResponse[t][m] * measuredCopy[m];
474
475 if (efficiency[t] > 0)
476 result[t] = value / efficiency[t];
477 else
478 result[t] = 0;
479 }
480
481 // draw intermediate result
482 /*
483 for (Int_t t=0; t<kMaxT; t++)
484 aResult->SetBinContent(t+1, result[t]);
485 aResult->SetMarkerStyle(20+i);
486 aResult->SetMarkerColor(2);
487 aResult->DrawCopy("P SAME HIST");
488 */
489
490 Double_t chi2LastIter = 0;
491 // regularization (simple smoothing)
492 for (Int_t t=kStartBin; t<kMaxT; t++)
493 {
494 Float_t newValue = 0;
495
496 // 0 bin excluded from smoothing
497 if (t > kStartBin+2 && t<kMaxT-1)
498 {
499 Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
500
501 // weight the average with the regularization parameter
502 newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
503 }
504 else
505 newValue = result[t];
506
507 // calculate chi2 (change from last iteration)
508 if (prior[t] > 1e-5)
509 {
510 Double_t diff = (prior[t] - newValue) / prior[t];
511 chi2LastIter += diff * diff;
512 }
513
514 prior[t] = newValue;
515 }
516 //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
517 //convergence->Fill(i+1, chi2LastIter);
518
519 if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit)
520 {
521 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);
522 break;
523 }
524 } // end of iterations
525
526 //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
527 //delete convergence;
528
529 for (Int_t t=0; t<kMaxT; t++)
530 aResult->SetBinContent(t+1, result[t]);
531
532 delete[] measuredCopy;
533 delete[] measuredError;
534 delete[] prior;
535 delete[] result;
536 delete[] efficiency;
537
538 for (Int_t i=0; i<kMaxT; i++)
539 {
540 delete[] response[i];
541 delete[] inverseResponse[i];
542 }
543 delete[] response;
544 delete[] inverseResponse;
545
546 return 0;
547
548 // ********
549 // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
550
551 /*printf("Calculating covariance matrix. This may take some time...\n");
552
553 // check if this is the right one...
554 TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
555
556 Int_t xBins = hInverseResponseBayes->GetNbinsX();
557 Int_t yBins = hInverseResponseBayes->GetNbinsY();
558
559 // calculate "unfolding matrix" Mij
560 Float_t matrixM[251][251];
561 for (Int_t i=1; i<=xBins; i++)
562 {
563 for (Int_t j=1; j<=yBins; j++)
564 {
565 if (fCurrentEfficiency->GetBinContent(i) > 0)
566 matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
567 else
568 matrixM[i-1][j-1] = 0;
569 }
570 }
571
572 Float_t* vectorn = new Float_t[yBins];
573 for (Int_t j=1; j<=yBins; j++)
574 vectorn[j-1] = fCurrentESD->GetBinContent(j);
575
576 // first part of covariance matrix, depends on input distribution n(E)
577 Float_t cov1[251][251];
578
579 Float_t nEvents = fCurrentESD->Integral(); // N
580
581 xBins = 20;
582 yBins = 20;
583
584 for (Int_t k=0; k<xBins; k++)
585 {
586 printf("In Cov1: %d\n", k);
587 for (Int_t l=0; l<yBins; l++)
588 {
589 cov1[k][l] = 0;
590
591 // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
592 for (Int_t j=0; j<yBins; j++)
593 cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
594 * (1.0 - vectorn[j] / nEvents);
595
596 // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
597 for (Int_t i=0; i<yBins; i++)
598 for (Int_t j=0; j<yBins; j++)
599 {
600 if (i == j)
601 continue;
602 cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
603 * vectorn[j] / nEvents;
604 }
605 }
606 }
607
608 printf("Cov1 finished\n");
609
610 TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
611 cov->Reset();
612
613 for (Int_t i=1; i<=xBins; i++)
614 for (Int_t j=1; j<=yBins; j++)
615 cov->SetBinContent(i, j, cov1[i-1][j-1]);
616
617 new TCanvas;
618 cov->Draw("COLZ");
619
620 // second part of covariance matrix, depends on response matrix
621 Float_t cov2[251][251];
622
623 // Cov[P(Er|Cu), P(Es|Cu)] term
624 Float_t covTerm[100][100][100];
625 for (Int_t r=0; r<yBins; r++)
626 for (Int_t u=0; u<xBins; u++)
627 for (Int_t s=0; s<yBins; s++)
628 {
629 if (r == s)
630 covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
631 * (1.0 - hResponse->GetBinContent(u+1, r+1));
632 else
633 covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
634 * hResponse->GetBinContent(u+1, s+1);
635 }
636
637 for (Int_t k=0; k<xBins; k++)
638 for (Int_t l=0; l<yBins; l++)
639 {
640 cov2[k][l] = 0;
641 printf("In Cov2: %d %d\n", k, l);
642 for (Int_t i=0; i<yBins; i++)
643 for (Int_t j=0; j<yBins; j++)
644 {
645 //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
646 // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
647 Float_t tmpCov = 0;
648 for (Int_t r=0; r<yBins; r++)
649 for (Int_t u=0; u<xBins; u++)
650 for (Int_t s=0; s<yBins; s++)
651 {
652 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
653 || hResponse->GetBinContent(u+1, i+1) == 0)
654 continue;
655
656 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
657 * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
658 * covTerm[r][u][s];
659 }
660
661 cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
662 }
663 }
664
665 printf("Cov2 finished\n");
666
667 for (Int_t i=1; i<=xBins; i++)
668 for (Int_t j=1; j<=yBins; j++)
669 cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
670
671 new TCanvas;
672 cov->Draw("COLZ");*/
673}
674
675//____________________________________________________________________
676Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
677{
678 // homogenity term for minuit fitting
679 // pure function of the parameters
680 // prefers constant function (pol0)
681
682 Double_t chi2 = 0;
683
684 for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
685 {
686 Double_t right = params[i];
687 Double_t left = params[i-1];
688
689 if (right != 0)
690 {
691 Double_t diff = 1 - left / right;
692 chi2 += diff * diff;
693 }
694 }
695
696 return chi2 / 100.0;
697}
698
699//____________________________________________________________________
700Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
701{
702 // homogenity term for minuit fitting
703 // pure function of the parameters
704 // prefers linear function (pol1)
705
706 Double_t chi2 = 0;
707
708 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
709 {
710 if (params[i-1] == 0)
711 continue;
712
713 Double_t right = params[i];
714 Double_t middle = params[i-1];
715 Double_t left = params[i-2];
716
717 Double_t der1 = (right - middle);
718 Double_t der2 = (middle - left);
719
720 Double_t diff = (der1 - der2) / middle;
721
722 chi2 += diff * diff;
723 }
724
725 return chi2;
726}
727
728//____________________________________________________________________
729Double_t AliUnfolding::RegularizationLog(TVectorD& params)
730{
731 // homogenity term for minuit fitting
732 // pure function of the parameters
733 // prefers linear function (pol1)
734
735 Double_t chi2 = 0;
736
737 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
738 {
739 if (params[i-1] == 0)
740 continue;
741
742 Double_t right = log(params[i]);
743 Double_t middle = log(params[i-1]);
744 Double_t left = log(params[i-2]);
745
746 Double_t der1 = (right - middle);
747 Double_t der2 = (middle - left);
748
749 Double_t diff = (der1 - der2) / middle;
750
751 chi2 += diff * diff;
752 }
753
754 return chi2 * 100;
755}
756
757//____________________________________________________________________
758Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
759{
760 // homogenity term for minuit fitting
761 // pure function of the parameters
762 // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
763 // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
764
765 Double_t chi2 = 0;
766
767 for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
768 {
769 Double_t right = params[i];
770 Double_t middle = params[i-1];
771 Double_t left = params[i-2];
772
773 Double_t der1 = (right - middle);
774 Double_t der2 = (middle - left);
775
776 Double_t diff = (der1 - der2);
777
778 chi2 += diff * diff;
779 }
780
781 return chi2 * 1e4;
782}
783
784//____________________________________________________________________
785Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
786{
787 // homogenity term for minuit fitting
788 // pure function of the parameters
789 // calculates entropy, from
790 // The method of reduced cross-entropy (M. Schmelling 1993)
791
792 Double_t paramSum = 0;
793
794 for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
795 paramSum += params[i];
796
797 Double_t chi2 = 0;
798 for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
799 {
800 //Double_t tmp = params[i] / paramSum;
801 Double_t tmp = params[i];
802 if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
803 {
804 chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
805 }
806 }
807
808 return 100.0 + chi2;
809}
810
811//____________________________________________________________________
812void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
813{
814 //
815 // fit function for minuit
816 // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
817 //
818
819 // d
820 static TVectorD paramsVector(fgMaxParams);
821 for (Int_t i=0; i<fgMaxParams; ++i)
822 paramsVector[i] = params[i] * params[i];
823
824 // calculate penalty factor
825 Double_t penaltyVal = 0;
826 switch (fgRegularizationType)
827 {
828 case kNone: break;
829 case kPol0: penaltyVal = RegularizationPol0(paramsVector); break;
830 case kPol1: penaltyVal = RegularizationPol1(paramsVector); break;
831 case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
832 case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
833 case kLog: penaltyVal = RegularizationLog(paramsVector); break;
834 }
835
836 //if (penaltyVal > 1e10)
837 // paramsVector2.Print();
838
839 penaltyVal *= fgRegularizationWeight;
840
841 // Ad
842 TVectorD measGuessVector(fgMaxInput);
843 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
844
845 // Ad - m
846 measGuessVector -= (*fgCurrentESDVector);
847
848 //measGuessVector.Print();
849
850 TVectorD copy(measGuessVector);
851
852 // (Ad - m) W
853 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
854 // normal way is like this:
855 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
856 // optimized way like this:
857 for (Int_t i=0; i<fgMaxInput; ++i)
858 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
859
860 //measGuessVector.Print();
861
862 // (Ad - m) W (Ad - m)
863 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
864 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
865 Double_t chi2FromFit = measGuessVector * copy * 1e6;
866
867 chi2 = chi2FromFit + penaltyVal;
868
869 static Int_t callCount = 0;
870 if ((callCount++ % 10000) == 0)
871 Printf("AliUnfolding::Chi2Function: Iteration %d: %f %f --> %f", callCount, chi2FromFit, penaltyVal, chi2);
872}
873
874//____________________________________________________________________
875void AliUnfolding::DrawGuess(Double_t *params)
876{
877 //
878 // draws residuals of solution suggested by params and effect of regularization
879 //
880
881 // d
882 static TVectorD paramsVector(fgMaxParams);
883 for (Int_t i=0; i<fgMaxParams; ++i)
884 paramsVector[i] = params[i] * params[i];
885
886 // Ad
887 TVectorD measGuessVector(fgMaxInput);
888 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
889
890 // Ad - m
891 measGuessVector -= (*fgCurrentESDVector);
892
893 TVectorD copy(measGuessVector);
894 //copy.Print();
895
896 // (Ad - m) W
897 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
898 // normal way is like this:
899 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
900 // optimized way like this:
901 for (Int_t i=0; i<fgMaxInput; ++i)
902 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
903
904 // (Ad - m) W (Ad - m)
905 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
906 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
907 //Double_t chi2FromFit = measGuessVector * copy * 1e6;
908
909 // draw residuals
910 TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
911 for (Int_t i=0; i<fgMaxInput; ++i)
912 residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
913 new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
914
915 // draw penalty
916 if (fgRegularizationType != kPol1) {
917 Printf("Drawing not supported");
918 return;
919 }
920
921 TH1* penalty = new TH1F("penalty", "penalty", fgMaxParams, -0.5, fgMaxParams - 0.5);
922
923 for (Int_t i=2+1; i<fgMaxParams; ++i)
924 {
925 if (params[i-1] == 0)
926 continue;
927
928 Double_t right = params[i];
929 Double_t middle = params[i-1];
930 Double_t left = params[i-2];
931
932 Double_t der1 = (right - middle);
933 Double_t der2 = (middle - left);
934
935 Double_t diff = (der1 - der2) / middle;
936
937 penalty->SetBinContent(i-1, diff * diff);
938
939 //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
940 }
941 new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
942}
943
944//____________________________________________________________________
945void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
946{
947 //
948 // fit function for minuit
949 // uses the TF1 stored in fgFitFunction
950 //
951
952 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
953 fgFitFunction->SetParameter(i, params[i]);
954
955 Double_t* params2 = new Double_t[fgMaxParams];
956
957 for (Int_t i=0; i<fgMaxParams; ++i)
958 params2[i] = fgFitFunction->Eval(i);
959
960 Chi2Function(unused1, unused2, chi2, params2, unused3);
961
962 delete[] params2;
963
964 if (fgDebug)
965 Printf("%f", chi2);
966}
967
968//____________________________________________________________________
969Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult)
970{
971 //
972 // correct spectrum using minuit chi2 method applying a functional fit
973 //
974
975 if (!fgFitFunction)
976 {
977 Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting.");
978 return -1;
979 }
980
981 SetChi2Regularization(kNone, 0);
982
983 SetStaticVariables(correlation, measured);
984
985 // Initialize TMinuit via generic fitter interface
986 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
987
988 minuit->SetFCN(TF1Function);
989 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
990 {
991 Double_t lower, upper;
992 fgFitFunction->GetParLimits(i, lower, upper);
993 minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper);
994 }
995
996 Double_t arglist[100];
997 arglist[0] = 0;
998 minuit->ExecuteCommand("SET PRINT", arglist, 1);
999 minuit->ExecuteCommand("MIGRAD", arglist, 0);
1000 //minuit->ExecuteCommand("MINOS", arglist, 0);
1001
1002 for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1003 fgFitFunction->SetParameter(i, minuit->GetParameter(i));
1004
1005 for (Int_t i=0; i<fgMaxParams; ++i)
1006 {
1007 Double_t value = fgFitFunction->Eval(i);
1008 if (fgDebug)
1009 Printf("%d : %f", i, value);
1010
1011 if (efficiency)
1012 {
1013 if (efficiency->GetBinContent(i+1) > 0)
1014 {
1015 value /= efficiency->GetBinContent(i+1);
1016 }
1017 else
1018 value = 0;
1019 }
1020 aResult->SetBinContent(i+1, value);
1021 aResult->SetBinError(i+1, 0);
1022 }
1023
1024 return 0;
1025}
1026
1027//____________________________________________________________________
1028void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
1029{
1030 // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins
1031 // The same limit is applied to the correlation
1032
1033 Int_t lastBin = 0;
1034 for (Int_t i=1; i<=measured->GetNbinsX(); ++i)
1035 {
1036 if (measured->GetBinContent(i) <= fgOverflowBinLimit)
1037 {
1038 lastBin = i;
1039 break;
1040 }
1041 }
1042
1043 if (lastBin == 0)
1044 {
1045 Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit);
1046 return;
1047 }
1048
1049 Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin);
1050
1051 TCanvas* canvas = 0;
1052
1053 if (fgDebug)
1054 {
1055 canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
1056 canvas->Divide(2, 2);
1057
1058 canvas->cd(1);
1059 measured->SetStats(kFALSE);
1060 measured->DrawCopy();
1061 gPad->SetLogy();
1062
1063 canvas->cd(2);
1064 correlation->SetStats(kFALSE);
1065 correlation->DrawCopy("COLZ");
1066 }
1067
1068 measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX()));
1069 for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i)
1070 {
1071 measured->SetBinContent(i, 0);
1072 measured->SetBinError(i, 0);
1073 }
1074 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1075 measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin)));
1076
1077 Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin));
1078
1079 for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1080 {
1081 correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY()));
1082 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1083 correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin)));
1084
1085 for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j)
1086 {
1087 correlation->SetBinContent(i, j, 0);
1088 correlation->SetBinError(i, j, 0);
1089 }
1090 }
1091
1092 Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!");
1093
1094 if (canvas)
1095 {
1096 canvas->cd(3);
1097 measured->DrawCopy();
1098 gPad->SetLogy();
1099
1100 canvas->cd(4);
1101 correlation->DrawCopy("COLZ");
1102 }
1103}