refactoring multiplicity analysis
[u/mrichter/AliRoot.git] / PWG0 / AliUnfolding.cxx
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
31 TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
32 TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
33 TVectorD* AliUnfolding::fgCurrentESDVector = 0;
34 TVectorD* AliUnfolding::fgEntropyAPriori = 0;
35
36 TF1* AliUnfolding::fgFitFunction = 0;
37
38 AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid;
39 Int_t AliUnfolding::fgMaxInput  = -1;  // bins in measured histogram
40 Int_t AliUnfolding::fgMaxParams = -1;  // bins in unfolded histogram = number of fit params
41 Float_t AliUnfolding::fgOverflowBinLimit = -1;
42
43 AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1;
44 Float_t AliUnfolding::fgRegularizationWeight = 10000;
45 Int_t AliUnfolding::fgSkipBinsBegin = 0;
46 Float_t AliUnfolding::fgMinuitStepSize = 0.1;                 // (usually not needed to be changed) step size in minimization
47 Bool_t AliUnfolding::fgNormalizeInput = kFALSE;                  // normalize input spectrum
48
49 Float_t AliUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
50 Int_t   AliUnfolding::fgBayesianIterations = 10;         // number of iterations in Bayesian method
51
52 Bool_t AliUnfolding::fgDebug = kFALSE;
53
54 ClassImp(AliUnfolding)
55
56 //____________________________________________________________________
57 void 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 //____________________________________________________________________
74 void 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 //____________________________________________________________________
84 void 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 //____________________________________________________________________
94 void 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 //____________________________________________________________________
104 void 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 //____________________________________________________________________
117 void 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 //____________________________________________________________________
130 void 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 //____________________________________________________________________
140 Int_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 //____________________________________________________________________
184 void 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 //____________________________________________________________________
247 Int_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 //____________________________________________________________________
359 Int_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   {
410     if (efficiency)
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 //____________________________________________________________________
676 Double_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 //____________________________________________________________________
700 Double_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 //____________________________________________________________________
729 Double_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 //____________________________________________________________________
758 Double_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 //____________________________________________________________________
785 Double_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 //____________________________________________________________________
812 void 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 //____________________________________________________________________
875 void 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 //____________________________________________________________________
945 void 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 //____________________________________________________________________
969 Int_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 //____________________________________________________________________
1028 void 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 }