Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / ANALYSIS / 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 #include <TExec.h>
31 #include "Riostream.h"
32 #include "TROOT.h"
33
34 using namespace std; //required for resolving the 'cout' symbol
35
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;
42
43 TAxis* AliUnfolding::fgUnfoldedAxis = 0;
44 TAxis* AliUnfolding::fgMeasuredAxis = 0;
45
46 TF1* AliUnfolding::fgFitFunction = 0;
47
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;
52
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 = 1000000;            // 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;
65
66 Float_t AliUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
67 Int_t   AliUnfolding::fgBayesianIterations = 10;          // number of iterations in Bayesian method
68
69 Bool_t AliUnfolding::fgDebug = kFALSE;
70
71 Int_t AliUnfolding::fgCallCount = 0;
72
73 Int_t AliUnfolding::fgPowern = 5;
74
75 Double_t AliUnfolding::fChi2FromFit = 0.;
76 Double_t AliUnfolding::fPenaltyVal  = 0.;
77 Double_t AliUnfolding::fAvgResidual = 0.;
78
79 Int_t AliUnfolding::fgPrintChi2Details = 0;
80
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;     
86
87 ClassImp(AliUnfolding)
88
89 //____________________________________________________________________
90 void AliUnfolding::SetUnfoldingMethod(MethodType methodType)
91 {
92   // set unfolding method
93   fgMethodType = methodType; 
94   
95   const char* name = 0;
96   switch (methodType)
97   {
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;
102   }
103   Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name);
104 }
105
106 //____________________________________________________________________
107 void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit) 
108
109   // enable the creation of a overflow bin that includes all statistics below the given limit
110   
111   fgOverflowBinLimit = overflowBinLimit; 
112   
113   Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit);
114 }
115
116 //____________________________________________________________________
117 void AliUnfolding::SetSkipBinsBegin(Int_t nBins)
118 {
119   // set number of skipped bins in regularization
120   
121   fgSkipBinsBegin = nBins;
122   
123   Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin);
124 }
125
126 //____________________________________________________________________
127 void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded) 
128
129   // set number of bins in the input (measured) distribution and in the unfolded distribution
130   fgMaxInput = nMeasured; 
131   fgMaxParams = nUnfolded; 
132   
133   if (fgCorrelationMatrix)
134   {
135     delete fgCorrelationMatrix;
136     fgCorrelationMatrix = 0;
137   }
138   if (fgCorrelationMatrixSquared)
139   {
140     fgCorrelationMatrixSquared = 0;
141     delete fgCorrelationMatrixSquared;
142   }
143   if (fgCorrelationCovarianceMatrix)
144   {
145     delete fgCorrelationCovarianceMatrix;
146     fgCorrelationCovarianceMatrix = 0;
147   }
148   if (fgCurrentESDVector)
149   {
150     delete fgCurrentESDVector;
151     fgCurrentESDVector = 0;
152   }
153   if (fgEntropyAPriori)
154   {
155     delete fgEntropyAPriori;
156     fgEntropyAPriori = 0;
157   }
158   if (fgEfficiency)
159   {
160     delete fgEfficiency;
161     fgEfficiency = 0;
162   }
163   if (fgUnfoldedAxis)
164   {
165     delete fgUnfoldedAxis;
166     fgUnfoldedAxis = 0;
167   }
168   if (fgMeasuredAxis)
169   {
170     delete fgMeasuredAxis;
171     fgMeasuredAxis = 0;
172   }
173   
174   Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
175 }
176
177 //____________________________________________________________________
178 void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight)
179 {
180   //
181   // sets the parameters for chi2 minimization
182   //
183
184   fgRegularizationType = type;
185   fgRegularizationWeight = weight;
186
187   Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight);
188 }
189
190 //____________________________________________________________________
191 void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
192 {
193   //
194   // sets the parameters for Bayesian unfolding
195   //
196
197   fgBayesianSmoothing = smoothing;
198   fgBayesianIterations = nIterations;
199
200   Printf("AliUnfolding::SetBayesianParameters --> Parameters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
201 }
202
203 //____________________________________________________________________
204 void AliUnfolding::SetFunction(TF1* function)
205 {
206   // set function for unfolding with a fit function
207   
208   fgFitFunction = function;
209   
210   Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar());
211 }
212
213 //____________________________________________________________________
214 Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
215 {
216   // unfolds with unfolding method fgMethodType
217   //
218   // parameters:
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
225   //
226   //  return code: see UnfoldWithMinuit/UnfoldWithBayesian/UnfoldWithFunction
227
228   if (fgMaxInput == -1)
229   {
230     Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
231     fgMaxInput = measured->GetNbinsX();
232   }
233   if (fgMaxParams == -1)
234   {
235     Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
236     fgMaxParams = measured->GetNbinsX();
237   }
238
239   if (fgOverflowBinLimit > 0)
240     CreateOverflowBin(correlation, measured);
241     
242   switch (fgMethodType)
243   {
244     case kInvalid:
245     {
246       Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting...");
247       return -1;
248     }
249     case kChi2Minimization:
250       return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check);
251     case kBayesian:
252       return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result);
253     case kFunction:
254       return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result);
255   }
256
257
258
259   return -1;
260 }
261
262 //____________________________________________________________________
263 void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency)
264 {
265   // fill static variables needed for minuit fit
266
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);
277   if (!fgEfficiency)
278     fgEfficiency = new TVectorD(fgMaxParams);
279   if (fgUnfoldedAxis)
280   {
281     delete fgUnfoldedAxis;
282     fgUnfoldedAxis = 0;
283   }
284   if (!fgUnfoldedAxis) 
285     fgUnfoldedAxis = new TAxis(*(correlation->GetXaxis()));
286   if (fgMeasuredAxis)
287   {
288     delete fgMeasuredAxis;
289     fgMeasuredAxis = 0;
290   }
291   if (!fgMeasuredAxis) 
292     fgMeasuredAxis = new TAxis(*(correlation->GetYaxis()));
293
294   fgCorrelationMatrix->Zero();
295   fgCorrelationCovarianceMatrix->Zero();
296   fgCurrentESDVector->Zero();
297   fgEntropyAPriori->Zero();
298
299   // normalize correction for given nPart
300   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
301   {
302     Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
303     if (sum <= 0)
304       continue;
305     Float_t maxValue = 0;
306     Int_t maxBin = -1;
307     for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
308     {
309       // find most probably value
310       if (maxValue < correlation->GetBinContent(i, j))
311       {
312         maxValue = correlation->GetBinContent(i, j);
313         maxBin = j;
314       }
315
316       // npart sum to 1
317       correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i));
318       correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
319
320       if (i <= fgMaxParams && j <= fgMaxInput)
321       {
322         (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
323         (*fgCorrelationMatrixSquared)(j-1, i-1) = correlation->GetBinContent(i, j) * correlation->GetBinContent(i, j);
324       }
325     }
326
327     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
328   }
329     
330   //normalize measured
331   Float_t smallestError = 1;
332   if (fgNormalizeInput)
333   {
334     Float_t sumMeasured = measured->Integral();
335     measured->Scale(1.0 / sumMeasured);
336     smallestError /= sumMeasured;
337   }
338   
339   for (Int_t i=0; i<fgMaxInput; ++i)
340   {
341     (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
342     if (measured->GetBinError(i+1) > 0)
343     {
344       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
345     }
346     else // in this case put error of 1, otherwise 0 bins are not added to the chi2...
347       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / smallestError / smallestError;
348
349     if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
350       (*fgCorrelationCovarianceMatrix)(i, i) = 0;
351     //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
352   }
353
354   // efficiency is expected to match bin width of result
355   for (Int_t i=0; i<fgMaxParams; ++i)
356   {
357     (*fgEfficiency)(i) = efficiency->GetBinContent(i+1);
358   }
359
360   if (correlation->GetNbinsX() != fgMaxParams || correlation->GetNbinsY() != fgMaxInput)
361     cout << "Response histo has incorrect dimensions; expect (" << fgMaxParams << ", " << fgMaxInput << "), got (" << correlation->GetNbinsX() << ", " << correlation->GetNbinsY() << ")" << endl;
362
363 }
364
365 //____________________________________________________________________
366 Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
367 {
368   //
369   // implementation of unfolding (internal function)
370   //
371   // unfolds <measured> using response from <correlation> and effiency <efficiency>
372   // output is in <result>
373   // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
374   //   negative values in initialConditions mean that the given parameter is fixed to the absolute of the value
375   // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
376   //
377   // returns minuit status (0 = success), (-1 when check was set)
378   //
379
380   SetStaticVariables(correlation, measured, efficiency);
381   
382   // Initialize TMinuit via generic fitter interface
383   Int_t params = fgMaxParams;
384   if (fgNotFoundEvents > 0)
385     params++;
386     
387   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, params);
388   Double_t arglist[100];
389   //  minuit->SetDefaultFitter("Minuit2");
390
391   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
392   arglist[0] = 0;
393   minuit->ExecuteCommand("SET PRINT", arglist, 1);
394
395   // however, enable warnings
396   //minuit->ExecuteCommand("SET WAR", arglist, 0);
397
398   // set minimization function
399   minuit->SetFCN(Chi2Function);
400
401   // set precision
402   minuit->SetPrecision(fgMinuitPrecision);
403
404   minuit->SetMaxIterations(fgMinuitMaxIterations);
405
406   minuit->ExecuteCommand("SET STRATEGY",&fgMinuitStrategy,1);
407
408   for (Int_t i=0; i<fgMaxParams; i++)
409     (*fgEntropyAPriori)[i] = 1;
410
411   // set initial conditions as a-priori distribution for MRX regularization
412   /*
413   for (Int_t i=0; i<fgMaxParams; i++)
414     if (initialConditions && initialConditions->GetBinContent(i+1) > 0)
415       (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
416   */
417
418   if (!initialConditions) {
419     initialConditions = measured;
420   } else {
421     Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions...");
422     //new TCanvas; initialConditions->DrawCopy();
423     if (fgNormalizeInput)
424       initialConditions->Scale(1.0 / initialConditions->Integral());
425   }
426
427   // extract minimum value from initial conditions (if we set a value to 0 it will stay 0)
428   Float_t minValue = 1e35;
429   if (fgMinimumInitialValueFix < 0)
430   {
431     for (Int_t i=0; i<fgMaxParams; ++i)
432     {
433       Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
434       if (initialConditions->GetBinContent(bin) > 0)
435         minValue = TMath::Min(minValue, (Float_t) initialConditions->GetBinContent(bin));
436     }
437   }
438   else
439     minValue = fgMinimumInitialValueFix;
440   
441   Double_t* results = new Double_t[fgMaxParams+1];
442   for (Int_t i=0; i<fgMaxParams; ++i)
443   {
444     Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
445     results[i] = initialConditions->GetBinContent(bin);
446
447     Bool_t fix = kFALSE;
448     if (results[i] < 0)
449     {
450       fix = kTRUE;
451       results[i] = -results[i];
452     }
453  
454     if (!fix && fgMinimumInitialValue && results[i] < minValue)
455       results[i] = minValue;
456       
457     // minuit sees squared values to prevent it from going negative...
458     results[i] = TMath::Sqrt(results[i]);
459
460     minuit->SetParameter(i, Form("param%d", i), results[i], (fix) ? 0 : fgMinuitStepSize, 0, 0);
461   }
462   if (fgNotFoundEvents > 0)
463   {
464     results[fgMaxParams] = efficiency->GetBinContent(1);
465     minuit->SetParameter(fgMaxParams, "vtx0", results[fgMaxParams], fgMinuitStepSize / 100, 0.01, 0.80);
466   }
467   
468   Int_t dummy = 0;
469   Double_t chi2 = 0;
470   Chi2Function(dummy, 0, chi2, results, 0);
471   printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2);
472
473   if (check)
474   {
475     DrawGuess(results);
476     delete[] results;
477     return -1;
478   }
479
480   // first param is number of iterations, second is precision....
481   arglist[0] = (float)fgMinuitMaxIterations;
482   //  arglist[1] = 1e-5;
483   //  minuit->ExecuteCommand("SET PRINT", arglist, 3);
484   //  minuit->ExecuteCommand("SCAN", arglist, 0);
485   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
486   Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status);
487   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
488   //minuit->ExecuteCommand("MINOS", arglist, 0);
489
490   if (fgNotFoundEvents > 0)
491   {
492     results[fgMaxParams] = minuit->GetParameter(fgMaxParams);
493     Printf("Efficiency for bin 0 changed from %f to %f", efficiency->GetBinContent(1), results[fgMaxParams]);
494     efficiency->SetBinContent(1, results[fgMaxParams]);
495   }
496   
497   for (Int_t i=0; i<fgMaxParams; ++i)
498   {
499     results[i] = minuit->GetParameter(i);
500     Double_t value = results[i] * results[i];
501    // error is : 2 * (relError on results[i]) * (value) = 2 * (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
502     Double_t error = 0;
503     if (TMath::IsNaN(minuit->GetParError(i)))
504       Printf("WARNING: Parameter %d error is nan", i);
505     else 
506       error = 2 * minuit->GetParError(i) * results[i];
507     
508     if (efficiency)
509     {   
510       //printf("value before efficiency correction: %f\n",value);
511       if (efficiency->GetBinContent(i+1) > 0)
512       {
513         value /= efficiency->GetBinContent(i+1);
514         error /= efficiency->GetBinContent(i+1);
515       }
516       else
517       {
518         value = 0;
519         error = 0;
520       }
521     }
522     //printf("value after efficiency correction: %f +/- %f\n",value,error);
523     result->SetBinContent(i+1, value);
524     result->SetBinError(i+1, error);
525   }
526
527   Int_t tmpCallCount = fgCallCount;
528   fgCallCount = 0; // needs to be 0 so that the Chi2Function prints its output
529   Chi2Function(dummy, 0, chi2, results, 0);
530   
531   Printf("AliUnfolding::UnfoldWithMinuit: iterations %d. Chi2 of final parameters is = %f", tmpCallCount, chi2);
532   
533   delete[] results;
534
535   return status;
536 }
537
538 //____________________________________________________________________
539 Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult)
540 {
541   //
542   // unfolds a spectrum using the Bayesian method
543   //
544   
545   if (measured->Integral() <= 0)
546   {
547     Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
548     return -1;
549   }
550
551   const Int_t kStartBin = 0;
552
553   Int_t kMaxM = fgMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
554   Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
555
556   // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
557   const Double_t kConvergenceLimit = kMaxT * 1e-6;
558
559   // store information in arrays, to increase processing speed (~ factor 5)
560   Double_t* measuredCopy = new Double_t[kMaxM];
561   Double_t* measuredError = new Double_t[kMaxM];
562   Double_t* prior = new Double_t[kMaxT];
563   Double_t* result = new Double_t[kMaxT];
564   Double_t* efficiency = new Double_t[kMaxT];
565   Double_t* binWidths = new Double_t[kMaxT];
566   Double_t* normResponse = new Double_t[kMaxT]; 
567
568   Double_t** response = new Double_t*[kMaxT];
569   Double_t** inverseResponse = new Double_t*[kMaxT];
570   for (Int_t i=0; i<kMaxT; i++)
571   {
572     response[i] = new Double_t[kMaxM];
573     inverseResponse[i] = new Double_t[kMaxM];
574     normResponse[i] = 0;
575   }
576
577   // for normalization
578   Float_t measuredIntegral = measured->Integral();
579   for (Int_t m=0; m<kMaxM; m++)
580   {
581     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
582     measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
583
584     for (Int_t t=0; t<kMaxT; t++)
585     {
586       response[t][m] = correlation->GetBinContent(t+1, m+1);
587       inverseResponse[t][m] = 0;
588       normResponse[t] += correlation->GetBinContent(t+1, m+1);
589     }
590   }
591
592   for (Int_t t=0; t<kMaxT; t++)
593   {
594     if (aEfficiency)
595     {
596       efficiency[t] = aEfficiency->GetBinContent(t+1);
597     }
598     else
599       efficiency[t] = 1;
600       
601     prior[t] = measuredCopy[t];
602     result[t] = 0;
603     binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1);
604
605     for (Int_t m=0; m<kMaxM; m++) { // Normalise response matrix
606       if (normResponse[t] != 0) 
607         response[t][m] /= normResponse[t];
608       else
609         Printf("AliUnfolding::UnfoldWithBayesian: Empty row,column in response matrix, for truth bin %d",t);
610     }
611   }
612
613   // pick prior distribution
614   if (initialConditions)
615   {
616     printf("Using different starting conditions...\n");
617     // for normalization
618     Float_t inputDistIntegral = initialConditions->Integral();
619     for (Int_t i=0; i<kMaxT; i++)
620       prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
621   }
622
623   //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
624   
625   //new TCanvas;
626   // unfold...
627   for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations < 0; i++)
628   {
629     if (fgDebug)
630       Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
631
632     // calculate IR from Bayes theorem
633     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
634
635     Double_t chi2Measured = 0;
636     for (Int_t m=0; m<kMaxM; m++)
637     {
638       Float_t norm = 0;
639       for (Int_t t = kStartBin; t<kMaxT; t++)
640         norm += response[t][m] * prior[t] * efficiency[t];
641
642       // calc. chi2: (measured - response * prior) / error
643       if (measuredError[m] > 0)
644       {
645         Double_t value = (measuredCopy[m] - norm) / measuredError[m];
646         chi2Measured += value * value;
647       }
648
649       if (norm > 0)
650       {
651         for (Int_t t = kStartBin; t<kMaxT; t++)
652           inverseResponse[t][m] = response[t][m] * prior[t] * efficiency[t] / norm;
653       }
654       else
655       {
656         for (Int_t t = kStartBin; t<kMaxT; t++)
657           inverseResponse[t][m] = 0;
658       }
659     }
660     //Printf("chi2Measured of the last prior is %e", chi2Measured);
661
662     for (Int_t t = kStartBin; t<kMaxT; t++)
663     {
664       Float_t value = 0;
665       for (Int_t m=0; m<kMaxM; m++)
666         value += inverseResponse[t][m] * measuredCopy[m];
667
668       if (efficiency[t] > 0)
669         result[t] = value / efficiency[t];
670       else
671         result[t] = 0;
672     }
673     
674     /* 
675     // draw intermediate result
676     for (Int_t t=0; t<kMaxT; t++)
677     {
678       aResult->SetBinContent(t+1, result[t]);
679     }
680     aResult->SetMarkerStyle(24+i);
681     aResult->SetMarkerColor(2);
682     aResult->DrawCopy((i == 0) ? "P" : "PSAME");
683     */
684  
685     Double_t chi2LastIter = 0;
686     // regularization (simple smoothing)
687     for (Int_t t=kStartBin; t<kMaxT; t++)
688     {
689       Float_t newValue = 0;
690       
691       // 0 bin excluded from smoothing
692       if (t > kStartBin+2 && t<kMaxT-1)
693       {
694         Float_t average = (result[t-1] / binWidths[t-1] + result[t] / binWidths[t] + result[t+1] / binWidths[t+1]) / 3 * binWidths[t];
695
696         // weight the average with the regularization parameter
697         newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
698       }
699       else
700         newValue = result[t];
701
702       // calculate chi2 (change from last iteration)
703       if (prior[t] > 1e-5)
704       {
705         Double_t diff = (prior[t] - newValue) / prior[t];
706         chi2LastIter += diff * diff;
707       }
708
709       prior[t] = newValue;
710     }
711     //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
712     //convergence->Fill(i+1, chi2LastIter);
713
714     if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit)
715     {
716       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);
717       break;
718     }
719   } // end of iterations
720
721   //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
722   //delete convergence;
723
724   Float_t factor = 1;
725   if (!fgNormalizeInput)
726     factor = measuredIntegral;
727   for (Int_t t=0; t<kMaxT; t++)
728     aResult->SetBinContent(t+1, result[t] * factor);
729
730   delete[] measuredCopy;
731   delete[] measuredError;
732   delete[] prior;
733   delete[] result;
734   delete[] efficiency;
735   delete[] binWidths;
736   delete[] normResponse;
737
738   for (Int_t i=0; i<kMaxT; i++)
739   {
740     delete[] response[i];
741     delete[] inverseResponse[i];
742   }
743   delete[] response;
744   delete[] inverseResponse;
745   
746   return 0;
747
748   // ********
749   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
750
751   /*printf("Calculating covariance matrix. This may take some time...\n");
752
753   // check if this is the right one...
754   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
755
756   Int_t xBins = hInverseResponseBayes->GetNbinsX();
757   Int_t yBins = hInverseResponseBayes->GetNbinsY();
758
759   // calculate "unfolding matrix" Mij
760   Float_t matrixM[251][251];
761   for (Int_t i=1; i<=xBins; i++)
762   {
763     for (Int_t j=1; j<=yBins; j++)
764     {
765       if (fCurrentEfficiency->GetBinContent(i) > 0)
766         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
767       else
768         matrixM[i-1][j-1] = 0;
769     }
770   }
771
772   Float_t* vectorn = new Float_t[yBins];
773   for (Int_t j=1; j<=yBins; j++)
774     vectorn[j-1] = fCurrentESD->GetBinContent(j);
775
776   // first part of covariance matrix, depends on input distribution n(E)
777   Float_t cov1[251][251];
778
779   Float_t nEvents = fCurrentESD->Integral(); // N
780
781   xBins = 20;
782   yBins = 20;
783
784   for (Int_t k=0; k<xBins; k++)
785   {
786     printf("In Cov1: %d\n", k);
787     for (Int_t l=0; l<yBins; l++)
788     {
789       cov1[k][l] = 0;
790
791       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
792       for (Int_t j=0; j<yBins; j++)
793         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
794           * (1.0 - vectorn[j] / nEvents);
795
796       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
797       for (Int_t i=0; i<yBins; i++)
798         for (Int_t j=0; j<yBins; j++)
799         {
800           if (i == j)
801             continue;
802           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
803             * vectorn[j] / nEvents;
804          }
805     }
806   }
807
808   printf("Cov1 finished\n");
809
810   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
811   cov->Reset();
812
813   for (Int_t i=1; i<=xBins; i++)
814     for (Int_t j=1; j<=yBins; j++)
815       cov->SetBinContent(i, j, cov1[i-1][j-1]);
816
817   new TCanvas;
818   cov->Draw("COLZ");
819
820   // second part of covariance matrix, depends on response matrix
821   Float_t cov2[251][251];
822
823   // Cov[P(Er|Cu), P(Es|Cu)] term
824   Float_t covTerm[100][100][100];
825   for (Int_t r=0; r<yBins; r++)
826     for (Int_t u=0; u<xBins; u++)
827       for (Int_t s=0; s<yBins; s++)
828       {
829         if (r == s)
830           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
831             * (1.0 - hResponse->GetBinContent(u+1, r+1));
832         else
833           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
834             * hResponse->GetBinContent(u+1, s+1);
835       }
836
837   for (Int_t k=0; k<xBins; k++)
838     for (Int_t l=0; l<yBins; l++)
839     {
840       cov2[k][l] = 0;
841       printf("In Cov2: %d %d\n", k, l);
842       for (Int_t i=0; i<yBins; i++)
843         for (Int_t j=0; j<yBins; j++)
844         {
845           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
846           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
847           Float_t tmpCov = 0;
848           for (Int_t r=0; r<yBins; r++)
849             for (Int_t u=0; u<xBins; u++)
850               for (Int_t s=0; s<yBins; s++)
851               {
852                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
853                   || hResponse->GetBinContent(u+1, i+1) == 0)
854                   continue;
855
856                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
857                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
858                   * covTerm[r][u][s];
859               }
860
861           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
862         }
863     }
864
865   printf("Cov2 finished\n");
866
867   for (Int_t i=1; i<=xBins; i++)
868     for (Int_t j=1; j<=yBins; j++)
869       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
870
871   new TCanvas;
872   cov->Draw("COLZ");*/
873 }
874
875 //____________________________________________________________________
876 Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
877 {
878   // homogenity term for minuit fitting
879   // pure function of the parameters
880   // prefers constant function (pol0)
881   //
882   // Does not take into account efficiency
883   Double_t chi2 = 0;
884
885   for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
886   {
887     Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
888     Double_t left   = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
889
890     if (left != 0)
891     {
892       Double_t diff = (right - left);
893       chi2 += diff * diff / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
894     }
895   }
896
897   return chi2 / 100.0;
898 }
899
900 //____________________________________________________________________
901 Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
902 {
903   // homogenity term for minuit fitting
904   // pure function of the parameters
905   // prefers linear function (pol1)
906   //
907   // Does not take into account efficiency
908   Double_t chi2 = 0;
909
910   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
911   {
912     if (params[i-1] == 0)
913       continue;
914
915     Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
916     Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
917     Double_t left   = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
918
919     Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
920     Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
921
922     //Double_t diff = (der1 - der2) / middle;
923     //chi2 += diff * diff;
924     chi2 += (der1 - der2) * (der1 - der2) / middle * fgUnfoldedAxis->GetBinWidth(i);
925   }
926
927   return chi2;
928 }
929
930 //____________________________________________________________________
931 Double_t AliUnfolding::RegularizationLog(TVectorD& params)
932 {
933   // homogenity term for minuit fitting
934   // pure function of the parameters
935   // prefers logarithmic function (log)
936   //
937   // Does not take into account efficiency
938
939   Double_t chi2 = 0;
940
941   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
942   {
943     if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
944      continue;
945
946     Double_t right  = log(params[i] / fgUnfoldedAxis->GetBinWidth(i+1));
947     Double_t middle = log(params[i-1] / fgUnfoldedAxis->GetBinWidth(i));
948     Double_t left   = log(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1));
949     
950     Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
951     Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
952     
953     //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
954
955     //if (fgCallCount == 0)
956     //  Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
957     chi2 += (der1 - der2) * (der1 - der2);// / error;
958   }
959
960   return chi2;
961 }
962
963 //____________________________________________________________________
964 Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
965 {
966   // homogenity term for minuit fitting
967   // pure function of the parameters
968   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
969   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
970   //
971   // Does not take into account efficiency
972
973   Double_t chi2 = 0;
974
975   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
976   {
977     Double_t right  = params[i];
978     Double_t middle = params[i-1];
979     Double_t left   = params[i-2];
980
981     Double_t der1 = (right - middle);
982     Double_t der2 = (middle - left);
983
984     Double_t diff = (der1 - der2);
985
986     chi2 += diff * diff;
987   }
988
989   return chi2 * 1e4;
990 }
991
992 //____________________________________________________________________
993 Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
994 {
995   // homogenity term for minuit fitting
996   // pure function of the parameters
997   // calculates entropy, from
998   // The method of reduced cross-entropy (M. Schmelling 1993)
999   //
1000   // Does not take into account efficiency
1001
1002   Double_t paramSum = 0;
1003   
1004   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
1005     paramSum += params[i];
1006
1007   Double_t chi2 = 0;
1008   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
1009   {
1010     Double_t tmp = params[i] / paramSum;
1011     //Double_t tmp = params[i];
1012     if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
1013     {
1014       chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
1015     }
1016     else
1017       chi2 += 100;
1018   }
1019
1020   return -chi2;
1021 }
1022
1023 //____________________________________________________________________
1024 Double_t AliUnfolding::RegularizationRatio(TVectorD& params)
1025 {
1026   // homogenity term for minuit fitting
1027   // pure function of the parameters
1028   //
1029   // Does not take into account efficiency
1030
1031   Double_t chi2 = 0;
1032
1033   for (Int_t i=5+fgSkipBinsBegin; i<fgMaxParams; ++i)
1034   {
1035     if (params[i-1] == 0 || params[i] == 0)
1036       continue;
1037
1038     Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1039     Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
1040     Double_t left   = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
1041     Double_t left2   = params[i-3] / fgUnfoldedAxis->GetBinWidth(i-2);
1042     Double_t left3   = params[i-4] / fgUnfoldedAxis->GetBinWidth(i-3);
1043     Double_t left4   = params[i-5] / fgUnfoldedAxis->GetBinWidth(i-4);
1044
1045     //Double_t diff = left / middle - middle / right;
1046     //Double_t diff = 2 * left / middle - middle / right - left2 / left;
1047     Double_t diff = 4 * left2 / left - middle / right - left / middle - left3 / left2 - left4 / left3;
1048     
1049     chi2 += diff * diff;// / middle;
1050   }
1051
1052   return chi2;
1053 }
1054
1055 //____________________________________________________________________
1056 Double_t AliUnfolding::RegularizationPowerLaw(TVectorD& params)
1057 {
1058   // homogenity term for minuit fitting
1059   // pure function of the parameters
1060   // prefers power law with n = -5
1061   //
1062   // Does not take into account efficiency
1063
1064   Double_t chi2 = 0;
1065
1066   Double_t right = 0.;
1067   Double_t middle = 0.;
1068   Double_t left = 0.;
1069
1070   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
1071   {
1072     if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1073       continue;
1074
1075     if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i-1) < 1e-8)
1076       continue;
1077     
1078     middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern);
1079
1080     if(middle>0) {
1081       right  = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i+1),fgPowern)/middle;
1082
1083       left   = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle;
1084       
1085       middle = 1.;
1086       
1087       Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1088       Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
1089     
1090       chi2 += (der1 - der2) * (der1 - der2)/ ( fgUnfoldedAxis->GetBinWidth(i+1)/2. + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)/2.)/( fgUnfoldedAxis->GetBinWidth(i+1)/2. + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)/2.);// / error;
1091       //   printf("i: %d chi2 = %f\n",i,chi2);
1092     }
1093
1094   }
1095
1096   return chi2;
1097 }
1098
1099 //____________________________________________________________________
1100 Double_t AliUnfolding::RegularizationLogLog(TVectorD& params)
1101 {
1102   // homogenity term for minuit fitting
1103   // pure function of the parameters
1104   // prefers a powerlaw (linear on a log-log scale)
1105   //
1106   // The calculation takes into account the efficiencies
1107
1108   Double_t chi2 = 0;
1109
1110   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
1111   {
1112     if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
1113      continue;
1114     if ((*fgEfficiency)(i-1) == 0 || (*fgEfficiency)(i) == 0 || (*fgEfficiency)(i-2) == 0)
1115      continue;
1116
1117     Double_t right  = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i+1));
1118     Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i));
1119     Double_t left   = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-1));
1120     
1121     Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) );
1122     Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) );
1123     
1124     double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.;
1125     Double_t dder = (der1-der2) / tmp;
1126
1127     chi2 += dder * dder;
1128   }
1129
1130   return chi2;
1131 }
1132
1133
1134
1135 //____________________________________________________________________
1136 void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
1137 {
1138   //
1139   // fit function for minuit
1140   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
1141   //
1142   
1143   // TODO use static members for the variables here to speed up processing (no construction/deconstruction)
1144
1145   // d = guess
1146   TVectorD paramsVector(fgMaxParams);
1147   for (Int_t i=0; i<fgMaxParams; ++i)
1148     paramsVector[i] = params[i] * params[i];
1149
1150   // calculate penalty factor
1151   Double_t penaltyVal = 0;
1152
1153   switch (fgRegularizationType)
1154   {
1155     case kNone:       break;
1156     case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
1157     case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
1158     case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
1159     case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
1160     case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
1161     case kRatio:      penaltyVal = RegularizationRatio(paramsVector); break;
1162     case kPowerLaw:   penaltyVal = RegularizationPowerLaw(paramsVector); break;
1163     case kLogLog:     penaltyVal = RegularizationLogLog(paramsVector); break;
1164   }
1165
1166   penaltyVal *= fgRegularizationWeight; //beta*PU
1167
1168   // Ad
1169   TVectorD measGuessVector(fgMaxInput);
1170   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1171
1172   // Ad - m
1173   measGuessVector -= (*fgCurrentESDVector);
1174   
1175 #if 0
1176   // new error calcuation using error on the guess
1177   
1178   // error from guess
1179   TVectorD errorGuessVector(fgMaxInput);
1180   //errorGuessVector = (*fgCorrelationMatrixSquared) * paramsVector;
1181   errorGuessVector = (*fgCorrelationMatrix) * paramsVector;
1182
1183   Double_t chi2FromFit = 0;
1184   for (Int_t i=0; i<fgMaxInput; ++i)
1185   {
1186     Float_t error = 1;
1187     if (errorGuessVector(i) > 0)
1188       error = errorGuessVector(i);
1189     chi2FromFit += measGuessVector(i) * measGuessVector(i) / error;
1190   }
1191
1192 #else
1193   // old error calcuation using the error on the measured
1194   TVectorD copy(measGuessVector);
1195
1196   // (Ad - m) W
1197   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1198   // normal way is like this:
1199   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1200   // optimized way like this:
1201   for (Int_t i=0; i<fgMaxInput; ++i)
1202     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1203
1204
1205   if (fgSkipBin0InChi2)
1206     measGuessVector[0] = 0;
1207
1208   // (Ad - m) W (Ad - m)
1209   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1210   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see UnfoldWithMinuit)
1211   Double_t chi2FromFit = measGuessVector * copy * 1e6;
1212 #endif
1213
1214   Double_t notFoundEventsConstraint = 0;
1215   Double_t currentNotFoundEvents = 0;
1216   Double_t errorNotFoundEvents = 0;
1217   
1218   if (fgNotFoundEvents > 0)
1219   {
1220     for (Int_t i=0; i<fgMaxParams; ++i)
1221     {
1222       Double_t eff = (1.0 / (*fgEfficiency)(i) - 1);
1223       if (i == 0)
1224         eff = (1.0 / params[fgMaxParams] - 1);
1225       currentNotFoundEvents += eff * paramsVector(i);
1226       errorNotFoundEvents += eff * eff * paramsVector(i); // error due to guess (paramsVector)
1227       if (i <= 3)
1228         errorNotFoundEvents += (eff * 0.03) * (eff * 0.03) * paramsVector(i) * paramsVector(i); // error on eff
1229       //      if ((fgCallCount % 10000) == 0)
1230         //Printf("%d %f %f %f", i, (*fgEfficiency)(i), paramsVector(i), currentNotFoundEvents);
1231     }
1232     errorNotFoundEvents += fgNotFoundEvents;
1233     // TODO add error on background, background estimate
1234     
1235     notFoundEventsConstraint = (currentNotFoundEvents - fgNotFoundEvents) * (currentNotFoundEvents - fgNotFoundEvents) / errorNotFoundEvents;
1236   }
1237   
1238   Float_t avg = 0;
1239   Float_t sum = 0;
1240   Float_t currentMult = 0;
1241   // try to match dn/deta
1242   for (Int_t i=0; i<fgMaxParams; ++i)
1243   {
1244     avg += paramsVector(i) * currentMult;
1245     sum += paramsVector(i);
1246     currentMult += fgUnfoldedAxis->GetBinWidth(i);
1247   }
1248   avg /= sum;
1249   Float_t chi2Avg = 0; //(avg - 3.73) * (avg - 3.73) * 100;
1250
1251   chi2 = chi2FromFit + penaltyVal + notFoundEventsConstraint + chi2Avg;
1252
1253   if ((fgCallCount++ % 1000) == 0)
1254   {
1255
1256     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);
1257
1258     //for (Int_t i=0; i<fgMaxInput; ++i)
1259     //  Printf("%d: %f", i, measGuessVector(i) * copy(i) * 1e6);
1260   }
1261
1262   fChi2FromFit = chi2FromFit;
1263   fPenaltyVal = penaltyVal;
1264 }
1265
1266 //____________________________________________________________________
1267 void AliUnfolding::MakePads() {
1268     TPad *presult = new TPad("presult","result",0,0.4,1,1);
1269     presult->SetNumber(1);
1270     presult->Draw();
1271     presult->SetLogy();
1272     TPad *pres = new TPad("pres","residuals",0,0.2,1,0.4);
1273     pres->SetNumber(2);
1274     pres->Draw();
1275     TPad *ppen = new TPad("ppen","penalty",0,0.0,1,0.2);
1276     ppen->SetNumber(3);
1277     ppen->Draw();
1278  
1279 }
1280 //____________________________________________________________________
1281 void AliUnfolding::DrawResults(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TCanvas *canv, Int_t reuseHists,TH1 *unfolded)
1282 {
1283   // Draw histograms of
1284   //   - Result folded with response
1285   //   - Penalty factors
1286   //   - Chisquare contributions
1287   // (Useful for debugging/sanity checks and the interactive unfolder)
1288   //
1289   // If a canvas pointer is given (canv != 0), it will be used for all
1290   // plots; 3 pads are made if needed.
1291
1292
1293   Int_t blankCanvas = 0;
1294   TVirtualPad *presult = 0;
1295   TVirtualPad *pres = 0; 
1296   TVirtualPad *ppen = 0;
1297   
1298   if (canv) {
1299     canv->cd();
1300     presult = canv->GetPad(1);
1301     pres = canv->GetPad(2);
1302     ppen = canv->GetPad(3); 
1303     if (presult == 0 || pres == 0 || ppen == 0) {
1304       canv->Clear();
1305       MakePads();
1306       blankCanvas = 1;
1307       presult = canv->GetPad(1);
1308       pres = canv->GetPad(2);
1309       ppen = canv->GetPad(3); 
1310     } 
1311   }
1312   else {
1313     presult = new TCanvas;
1314     pres = new TCanvas;
1315     ppen = new TCanvas;
1316   }
1317
1318
1319   if (fgMaxInput == -1)
1320   {
1321     Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
1322     fgMaxInput = measured->GetNbinsX();
1323   }
1324   if (fgMaxParams == -1)
1325   {
1326     Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
1327     fgMaxParams = initialConditions->GetNbinsX();
1328   }
1329
1330   if (fgOverflowBinLimit > 0)
1331     CreateOverflowBin(correlation, measured);
1332
1333   // Uses Minuit methods
1334
1335   SetStaticVariables(correlation, measured, efficiency);
1336
1337   Double_t* params = new Double_t[fgMaxParams+1];
1338   for (Int_t i=0; i<fgMaxParams; ++i)
1339   {
1340     params[i] = initialConditions->GetBinContent(i+1) * efficiency->GetBinContent(i+1);
1341
1342     Bool_t fix = kFALSE;
1343     if (params[i] < 0)
1344     {
1345       fix = kTRUE;
1346       params[i] = -params[i];
1347     }
1348     params[i]=TMath::Sqrt(params[i]);
1349   } 
1350
1351   Double_t chi2;
1352   Int_t dummy;
1353
1354   //fgPrintChi2Details = kTRUE;
1355   fgCallCount = 0; // To make sure that Chi2Function prints the components
1356   Chi2Function(dummy, 0, chi2, params, 0);
1357
1358   presult->cd();
1359   TH1 *meas2 = (TH1*)measured->Clone("meas2");
1360   meas2->SetTitle("meas2");
1361   meas2->SetName("meas2");
1362   meas2->SetLineColor(1);
1363   meas2->SetLineWidth(2);
1364   meas2->SetMarkerColor(meas2->GetLineColor());
1365   meas2->SetMarkerStyle(20);
1366   Float_t scale = unfolded->GetXaxis()->GetBinWidth(1)/meas2->GetXaxis()->GetBinWidth(1);
1367   meas2->Scale(scale);
1368   if (blankCanvas)
1369     meas2->DrawCopy();
1370   else
1371     meas2->DrawCopy("same");
1372   //meas2->GetXaxis()->SetLimits(0,fgMaxInput);
1373   meas2->SetBit(kCannotPick);
1374   DrawGuess(params, presult, pres, ppen, reuseHists,unfolded);
1375   delete [] params;
1376 }
1377 //____________________________________________________________________
1378 void AliUnfolding::RedrawInteractive() {
1379   // 
1380   // Helper function for interactive unfolding
1381   //
1382   DrawResults(fghCorrelation,fghEfficiency,fghMeasured,fghUnfolded,fgCanvas,1,fghUnfolded);
1383 }
1384 //____________________________________________________________________
1385 void AliUnfolding::InteractiveUnfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions) 
1386 {
1387   //
1388   // Function to do interactive unfolding
1389   // A canvas is drawn with the unfolding result
1390   // Change the histogram with your mouse and all histograms 
1391   // will be updated automatically
1392
1393   fgCanvas = new TCanvas("UnfoldingCanvas","Interactive unfolding",500,800);
1394   fgCanvas->cd();
1395
1396   MakePads();
1397
1398   if (fghUnfolded) {
1399     delete fghUnfolded; fghUnfolded = 0;
1400   }
1401   
1402   gROOT->SetEditHistograms(kTRUE);
1403
1404   fghUnfolded = new TH1F("AliUnfoldingfghUnfolded","Unfolded result (Interactive unfolder",efficiency->GetNbinsX(),efficiency->GetXaxis()->GetXmin(),efficiency->GetXaxis()->GetXmax());
1405
1406   fghCorrelation = correlation;
1407   fghEfficiency = efficiency;
1408   fghMeasured = measured;
1409
1410   if(fgMinuitMaxIterations>0)
1411     UnfoldWithMinuit(correlation,efficiency,measured,initialConditions,fghUnfolded,kFALSE);
1412   else
1413     fghUnfolded = initialConditions;
1414
1415   fghUnfolded->SetLineColor(2);
1416   fghUnfolded->SetMarkerColor(2);
1417   fghUnfolded->SetLineWidth(2);
1418
1419
1420   fgCanvas->cd(1);
1421   fghUnfolded->Draw();
1422   DrawResults(correlation,efficiency,measured,fghUnfolded,fgCanvas,kFALSE,fghUnfolded);
1423
1424   TExec *execRedraw = new TExec("redraw","AliUnfolding::RedrawInteractive()");
1425   fghUnfolded->GetListOfFunctions()->Add(execRedraw);
1426 }
1427 //____________________________________________________________________
1428 void AliUnfolding::DrawGuess(Double_t *params, TVirtualPad *pfolded, TVirtualPad *pres, TVirtualPad *ppen, Int_t reuseHists,TH1* unfolded)
1429 {
1430   //
1431   // draws residuals of solution suggested by params and effect of regularization
1432   //
1433
1434   if (pfolded == 0)
1435     pfolded = new TCanvas;
1436   if (ppen == 0)
1437     ppen = new TCanvas;
1438   if (pres == 0)
1439     pres = new TCanvas;
1440   
1441   // d
1442   TVectorD paramsVector(fgMaxParams);
1443   for (Int_t i=0; i<fgMaxParams; ++i)
1444     paramsVector[i] = params[i] * params[i];
1445
1446   // Ad
1447   TVectorD measGuessVector(fgMaxInput);
1448   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1449
1450   TH1* folded = 0;
1451   if (reuseHists) 
1452     folded = dynamic_cast<TH1*>(gROOT->FindObject("__hfolded"));
1453   if (!reuseHists || folded == 0) {
1454     if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins
1455       folded = new TH1F("__hfolded","Folded histo from AliUnfolding",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray());
1456     else
1457       folded = new TH1F("__hfolded","Folded histo from AliUnfolding",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(),fgMeasuredAxis->GetXmax());
1458   }
1459
1460   folded->SetBit(kCannotPick);
1461   folded->SetLineColor(4);
1462   folded->SetLineWidth(2);
1463
1464   for (Int_t ibin =0; ibin < fgMaxInput; ibin++) 
1465     folded->SetBinContent(ibin+1, measGuessVector[ibin]);
1466
1467   Float_t scale = unfolded->GetXaxis()->GetBinWidth(1)/folded->GetXaxis()->GetBinWidth(1);
1468   folded->Scale(scale);
1469
1470   pfolded->cd();
1471
1472   folded->Draw("same");
1473
1474   // Ad - m
1475   measGuessVector -= (*fgCurrentESDVector);
1476
1477   TVectorD copy(measGuessVector);
1478
1479   // (Ad - m) W
1480   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1481   // normal way is like this:
1482   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1483   // optimized way like this:
1484   for (Int_t i=0; i<fgMaxInput; ++i)
1485     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1486
1487   // (Ad - m) W (Ad - m)
1488   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1489   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
1490   //Double_t chi2FromFit = measGuessVector * copy * 1e6;
1491
1492   // draw residuals
1493   // Double_t pTarray[fgMaxParams+1];
1494   // for(int i=0; i<fgMaxParams; i++) {
1495   //   pTarray[i] = fgUnfoldedAxis->GetBinCenter(i)-0.5*fgUnfoldedAxis->GetBinWidth(i);
1496   // }
1497   //  pTarray[fgMaxParams] = fgUnfoldedAxis->GetBinCenter(fgMaxParams-1)+0.5*fgUnfoldedAxis->GetBinWidth(fgMaxParams-1);
1498   //  TH1* residuals = new TH1F("residuals", "residuals", fgMaxParams,pTarray);
1499   //  TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
1500   // for (Int_t i=0; i<fgMaxInput; ++i)
1501   //   residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);7
1502   TH1* residuals = GetResidualsPlot(params);
1503   residuals->GetXaxis()->SetTitleSize(0.06);
1504   residuals->GetXaxis()->SetTitleOffset(0.7);
1505   residuals->GetXaxis()->SetLabelSize(0.07);
1506   residuals->GetYaxis()->SetTitleSize(0.08);
1507   residuals->GetYaxis()->SetTitleOffset(0.5);
1508   residuals->GetYaxis()->SetLabelSize(0.06);
1509   pres->cd(); residuals->DrawCopy(); gPad->SetLogy();
1510     
1511
1512   // draw penalty
1513   TH1* penalty = GetPenaltyPlot(params);
1514   penalty->GetXaxis()->SetTitleSize(0.06);
1515   penalty->GetXaxis()->SetTitleOffset(0.7);
1516   penalty->GetXaxis()->SetLabelSize(0.07);
1517   penalty->GetYaxis()->SetTitleSize(0.08);
1518   penalty->GetYaxis()->SetTitleOffset(0.5);
1519   penalty->GetYaxis()->SetLabelSize(0.06);
1520   ppen->cd(); penalty->DrawCopy(); gPad->SetLogy();
1521 }
1522
1523 //____________________________________________________________________
1524 TH1* AliUnfolding::GetResidualsPlot(TH1* corrected)
1525 {
1526   //
1527   // MvL: THIS MUST BE INCORRECT. 
1528   // Need heff to calculate params from TH1 'corrected'
1529   //
1530
1531   //
1532   // fill residuals histogram of solution suggested by params and effect of regularization
1533   //
1534
1535   Double_t* params = new Double_t[fgMaxParams];
1536   for (Int_t i=0; i<fgMaxParams; i++)
1537     params[i] = 0;
1538
1539   for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
1540     params[i] = TMath::Sqrt(TMath::Abs(corrected->GetBinContent(i+1)*(*fgEfficiency)(i)));
1541
1542   TH1 * plot = GetResidualsPlot(params);
1543   delete [] params;
1544   return plot;
1545 }
1546
1547 //____________________________________________________________________
1548 TH1* AliUnfolding::GetResidualsPlot(Double_t* params)
1549 {
1550   //
1551   // fill residuals histogram of solution suggested by params and effect of regularization
1552   //
1553
1554   // d
1555   TVectorD paramsVector(fgMaxParams);
1556   for (Int_t i=0; i<fgMaxParams; ++i)
1557     paramsVector[i] = params[i] * params[i];
1558
1559   // Ad
1560   TVectorD measGuessVector(fgMaxInput);
1561   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1562
1563   // Ad - m
1564   measGuessVector -= (*fgCurrentESDVector);
1565
1566   TVectorD copy(measGuessVector);
1567
1568   // (Ad - m) W
1569   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1570   // normal way is like this:
1571   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1572   // optimized way like this:
1573   for (Int_t i=0; i<fgMaxInput; ++i)
1574     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1575
1576   // (Ad - m) W (Ad - m)
1577   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1578   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
1579   //Double_t chi2FromFit = measGuessVector * copy * 1e6;
1580
1581   // draw residuals
1582   TH1* residuals = 0;
1583   if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins
1584     residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray());
1585   else
1586     residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(), fgMeasuredAxis->GetXmax());
1587   //  TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
1588
1589   Double_t sumResiduals = 0.; 
1590   for (Int_t i=0; i<fgMaxInput; ++i) {
1591     residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
1592     sumResiduals += measGuessVector(i) * copy(i) * 1e6;
1593   }
1594   fAvgResidual = sumResiduals/(double)fgMaxInput;
1595  
1596   return residuals;
1597 }
1598
1599 //____________________________________________________________________
1600 TH1* AliUnfolding::GetPenaltyPlot(TH1* corrected)
1601 {
1602   // draws the penalty factors as function of multiplicity of the current selected regularization
1603
1604   Double_t* params = new Double_t[fgMaxParams];
1605   for (Int_t i=0; i<fgMaxParams; i++)
1606     params[i] = 0;
1607
1608   for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
1609     params[i] = (*fgEfficiency)(i)*corrected->GetBinContent(i+1);
1610   
1611   TH1* penalty = GetPenaltyPlot(params);
1612   
1613   delete[] params;
1614   
1615   return penalty;
1616 }
1617
1618 //____________________________________________________________________
1619 TH1* AliUnfolding::GetPenaltyPlot(Double_t* params)
1620 {
1621   // draws the penalty factors as function of multiplicity of the current selected regularization
1622   
1623   //TH1* penalty = new TH1F("penalty", ";unfolded multiplicity;penalty factor", fgMaxParams, -0.5, fgMaxParams - 0.5);
1624   //  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) );
1625
1626   TH1* penalty = 0;
1627   if (fgUnfoldedAxis->GetXbins()->GetArray())
1628     penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXbins()->GetArray());
1629   else
1630     penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXmin(),fgUnfoldedAxis->GetXmax());
1631
1632   for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
1633   {
1634     Double_t diff = 0;
1635     if (fgRegularizationType == kPol0)
1636     {
1637       Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1638       Double_t left   = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
1639
1640       if (left != 0)
1641       {
1642         Double_t diffTmp = (right - left);
1643         diff = diffTmp * diffTmp / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2) / 100;
1644       }
1645     } 
1646     if (fgRegularizationType == kPol1 && i > 1+fgSkipBinsBegin) 
1647     {
1648       if (params[i-1] == 0)
1649         continue;
1650
1651       Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1652       Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
1653       Double_t left   = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
1654
1655       Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1656       Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
1657
1658       diff = (der1 - der2) * (der1 - der2) / middle;
1659     }
1660
1661     if (fgRegularizationType == kLog && i > 1+fgSkipBinsBegin) 
1662     {
1663       if (params[i-1] == 0)
1664         continue;
1665   
1666       Double_t right  = log(params[i]);
1667       Double_t middle = log(params[i-1]);
1668       Double_t left   = log(params[i-2]);
1669       
1670       Double_t der1 = (right - middle);
1671       Double_t der2 = (middle - left);
1672   
1673       //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
1674       //Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
1675        
1676       diff = (der1 - der2) * (der1 - der2);// / error;
1677     }
1678     if (fgRegularizationType == kCurvature && i > 1+fgSkipBinsBegin)
1679     {
1680       Double_t right  = params[i];    // params are sqrt
1681       Double_t middle = params[i-1];
1682       Double_t left   = params[i-2];
1683       
1684       Double_t der1 = (right - middle)/0.5/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i));
1685       Double_t der2 = (middle - left)/0.5/(fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i+1));
1686       
1687       diff = (der1 - der2)/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1))*3.0;
1688       diff = 1e4*diff*diff;
1689     }
1690     if (fgRegularizationType == kPowerLaw && i > 1+fgSkipBinsBegin) 
1691     {
1692
1693       if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1694         continue;
1695       
1696       if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8)
1697         continue;
1698       
1699       double middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern);
1700
1701       if(middle>0) {
1702         double right  = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i+1),fgPowern)/middle;
1703         
1704         double left   = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle;
1705         
1706         middle = 1.;
1707         
1708         Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1709         Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
1710         
1711         diff = (der1 - der2) * (der1 - der2);// / error;
1712       }
1713     }
1714
1715     if (fgRegularizationType == kLogLog && i > 1+fgSkipBinsBegin) 
1716     {
1717
1718       if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1719         continue;
1720
1721       Double_t right  = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i+1));
1722       Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i));
1723       Double_t left   = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-1));
1724       
1725       Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) );
1726       Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) );
1727       
1728       double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.;
1729       Double_t dder = (der1-der2) / tmp;
1730       
1731       diff = dder * dder;
1732     }
1733     
1734     penalty->SetBinContent(i, diff*fgRegularizationWeight);
1735     
1736     //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
1737   }
1738   
1739   return penalty;
1740 }
1741
1742 //____________________________________________________________________
1743 void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
1744 {
1745   //
1746   // fit function for minuit
1747   // uses the TF1 stored in fgFitFunction
1748   //
1749
1750   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1751     fgFitFunction->SetParameter(i, params[i]);
1752
1753   Double_t* params2 = new Double_t[fgMaxParams];
1754
1755   for (Int_t i=0; i<fgMaxParams; ++i)
1756     params2[i] = fgFitFunction->Eval(i);
1757
1758   Chi2Function(unused1, unused2, chi2, params2, unused3);
1759   
1760   delete[] params2;
1761
1762   if (fgDebug)
1763     Printf("%f", chi2);
1764 }
1765
1766 //____________________________________________________________________
1767 Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult)
1768 {
1769   //
1770   // correct spectrum using minuit chi2 method applying a functional fit
1771   //
1772   
1773   if (!fgFitFunction)
1774   {
1775     Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting.");
1776     return -1;
1777   }    
1778   
1779   SetChi2Regularization(kNone, 0);
1780   
1781   SetStaticVariables(correlation, measured, efficiency);
1782   
1783   // Initialize TMinuit via generic fitter interface
1784   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
1785
1786   minuit->SetFCN(TF1Function);
1787   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1788   {
1789     Double_t lower, upper;
1790     fgFitFunction->GetParLimits(i, lower, upper);
1791     minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper);
1792   }
1793
1794   Double_t arglist[100];
1795   arglist[0] = 0;
1796   minuit->ExecuteCommand("SET PRINT", arglist, 1);
1797   minuit->ExecuteCommand("SCAN", arglist, 0);
1798   minuit->ExecuteCommand("MIGRAD", arglist, 0);
1799   //minuit->ExecuteCommand("MINOS", arglist, 0);
1800
1801   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1802     fgFitFunction->SetParameter(i, minuit->GetParameter(i));
1803
1804   for (Int_t i=0; i<fgMaxParams; ++i)
1805   {
1806     Double_t value = fgFitFunction->Eval(i);
1807     if (fgDebug)
1808       Printf("%d : %f", i, value);
1809     
1810     if (efficiency)
1811     {
1812       if (efficiency->GetBinContent(i+1) > 0)
1813       {
1814         value /= efficiency->GetBinContent(i+1);
1815       }
1816       else
1817         value = 0;
1818     }
1819     aResult->SetBinContent(i+1, value);
1820     aResult->SetBinError(i+1, 0);
1821   }
1822   
1823   return 0;
1824 }
1825
1826 //____________________________________________________________________
1827 void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
1828 {
1829   // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins
1830   // The same limit is applied to the correlation  
1831   
1832   Int_t lastBin = 0;
1833   for (Int_t i=1; i<=measured->GetNbinsX(); ++i)
1834   {
1835     if (measured->GetBinContent(i) <= fgOverflowBinLimit)
1836     {
1837       lastBin = i;
1838       break;
1839     }
1840   }
1841   
1842   if (lastBin == 0)
1843   {
1844     Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit);
1845     return;
1846   }
1847   
1848   Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin);
1849   
1850   TCanvas* canvas = 0;
1851
1852   if (fgDebug)
1853   {
1854     canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
1855     canvas->Divide(2, 2);
1856
1857     canvas->cd(1);
1858     measured->SetStats(kFALSE);
1859     measured->DrawCopy();
1860     gPad->SetLogy();
1861
1862     canvas->cd(2);
1863     correlation->SetStats(kFALSE);
1864     correlation->DrawCopy("COLZ");
1865   }
1866
1867   measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX()));
1868   for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i)
1869   {
1870     measured->SetBinContent(i, 0);
1871     measured->SetBinError(i, 0);
1872   }
1873   // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1874   measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin)));
1875
1876   Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin));
1877
1878   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1879   {
1880     correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY()));
1881     // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1882     correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin)));
1883
1884     for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j)
1885     {
1886       correlation->SetBinContent(i, j, 0);
1887       correlation->SetBinError(i, j, 0);
1888     }
1889   }
1890
1891   Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!");
1892
1893   if (canvas)
1894   {
1895     canvas->cd(3);
1896     measured->DrawCopy();
1897     gPad->SetLogy();
1898
1899     canvas->cd(4);
1900     correlation->DrawCopy("COLZ");
1901   }
1902 }
1903
1904 Int_t AliUnfolding::UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result)
1905 {
1906   // unfolds and assigns bias as errors with Eq. 19 of Cowan, "a survey of unfolding methods for particle physics"
1907   // b_i = sum_j dmu_i/dn_j (nu_j - n_j) with nu_j as folded guess, n_j as data
1908   // dmu_i/dn_j is found numerically by changing the bin content and re-unfolding
1909   //
1910   // return code: 0 (success) -1 (error: from Unfold(...) )
1911
1912   if (Unfold(correlation, efficiency, measured, initialConditions, result) != 0)
1913     return -1;
1914
1915   TMatrixD matrix(fgMaxInput, fgMaxParams);
1916   
1917   TH1* newResult[4];
1918   for (Int_t loop=0; loop<4; loop++)
1919     newResult[loop] = (TH1*) result->Clone(Form("newresult_%d", loop));
1920     
1921   // change bin-by-bin and built matrix of effects
1922   for (Int_t m=0; m<fgMaxInput; m++)
1923   {
1924     if (measured->GetBinContent(m+1) < 1)
1925       continue;
1926       
1927     for (Int_t loop=0; loop<4; loop++)
1928     {
1929       //Printf("%d %d", i, loop);
1930       Float_t factor = 1;
1931       switch (loop)
1932       {
1933         case 0: factor = 0.5; break;
1934         case 1: factor = -0.5; break;
1935         case 2: factor = 1; break;
1936         case 3: factor = -1; break;
1937         default: return -1;
1938       }
1939       
1940       TH1* measuredClone = (TH1*) measured->Clone("measuredClone");
1941   
1942       measuredClone->SetBinContent(m+1, measured->GetBinContent(m+1) + factor * TMath::Sqrt(measured->GetBinContent(m+1)));
1943       //new TCanvas; measuredClone->Draw("COLZ");
1944     
1945       newResult[loop]->Reset();
1946       Unfold(correlation, efficiency, measuredClone, measuredClone, newResult[loop]);
1947       // WARNING if we leave here, then nothing is calculated
1948         //return -1;
1949       
1950       delete measuredClone;
1951     }
1952     
1953     for (Int_t t=0; t<fgMaxParams; t++)
1954     {
1955       // one value
1956       //matrix(i, j) = (result->GetBinContent(j+1) - newResult->GetBinContent(j+1)) / TMath::Sqrt(changedHist->GetBinContent(1, i+1));
1957       
1958       // four values from the derivate (procedure taken from ROOT -- suggestion by Ruben)
1959       // = 1/2D * [ 8 (f(D/2) - f(-D/2)) - (f(D)-f(-D)) ]/3
1960       
1961       /*
1962       // check formula
1963       measured->SetBinContent(1, m+1, 100);
1964       newResult[0]->SetBinContent(t+1, 5 * 0.5 + 10);
1965       newResult[1]->SetBinContent(t+1, 5 * -0.5 + 10);
1966       newResult[2]->SetBinContent(t+1, 5 * 1 + 10);
1967       newResult[3]->SetBinContent(t+1, 5 * -1 + 10);
1968       */
1969       
1970       matrix(m, t) = 0.5 / TMath::Sqrt(measured->GetBinContent(m+1)) * 
1971         ( 8. * (newResult[0]->GetBinContent(t+1) - newResult[1]->GetBinContent(t+1)) - 
1972               (newResult[2]->GetBinContent(t+1) - newResult[3]->GetBinContent(t+1))
1973         ) / 3;
1974     }
1975   }
1976   
1977   for (Int_t loop=0; loop<4; loop++)
1978     delete newResult[loop];
1979   
1980   // ...calculate folded guess  
1981   TH1* convoluted = (TH1*) measured->Clone("convoluted");
1982   convoluted->Reset();
1983   for (Int_t m=0; m<fgMaxInput; m++)
1984   {
1985     Float_t value = 0;
1986     for (Int_t t = 0; t<fgMaxParams; t++)
1987     {
1988       Float_t tmp = correlation->GetBinContent(t+1, m+1) * result->GetBinContent(t+1);
1989       if (efficiency)
1990         tmp *= efficiency->GetBinContent(t+1);
1991       value += tmp;
1992     }
1993     convoluted->SetBinContent(m+1, value);
1994   }
1995   
1996   // rescale
1997   convoluted->Scale(measured->Integral() / convoluted->Integral());
1998   
1999   //new TCanvas; convoluted->Draw(); measured->Draw("SAME"); measured->SetLineColor(2);
2000   //return;
2001   
2002   // difference
2003   convoluted->Add(measured, -1);
2004   
2005   // ...and the bias
2006   for (Int_t t = 0; t<fgMaxParams; t++)
2007   {
2008     Double_t error = 0;
2009     for (Int_t m=0; m<fgMaxInput; m++)
2010       error += matrix(m, t) * convoluted->GetBinContent(m+1); 
2011     result->SetBinError(t+1, error);
2012   }
2013   
2014   //new TCanvas; result->Draw(); gPad->SetLogy();
2015   
2016   return 0;
2017 }