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