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