Updated AliUnfolding class. The regularization schemas have changed taking into accou...
[u/mrichter/AliRoot.git] / PWG0 / AliUnfolding.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: AliUnfolding.cxx 31168 2009-02-23 15:18:45Z jgrosseo $ */
17
18 // This class allows 1-dimensional unfolding.
19 // Methods that are implemented are chi2 minimization and bayesian unfolding.
20 //
21 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
22
23 #include "AliUnfolding.h"
24 #include <TH1F.h>
25 #include <TH2F.h>
26 #include <TVirtualFitter.h>
27 #include <TMath.h>
28 #include <TCanvas.h>
29 #include <TF1.h>
30
31 TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
32 TMatrixD* AliUnfolding::fgCorrelationMatrixSquared = 0;
33 TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
34 TVectorD* AliUnfolding::fgCurrentESDVector = 0;
35 TVectorD* AliUnfolding::fgEntropyAPriori = 0;
36 TVectorD* AliUnfolding::fgEfficiency = 0;
37 TVectorD* AliUnfolding::fgBinWidths = 0;
38
39 TF1* AliUnfolding::fgFitFunction = 0;
40
41 AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid;
42 Int_t AliUnfolding::fgMaxInput  = -1;  // bins in measured histogram
43 Int_t AliUnfolding::fgMaxParams = -1;  // bins in unfolded histogram = number of fit params
44 Float_t AliUnfolding::fgOverflowBinLimit = -1;
45
46 AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1;
47 Float_t AliUnfolding::fgRegularizationWeight = 10000;
48 Int_t AliUnfolding::fgSkipBinsBegin = 0;
49 Float_t AliUnfolding::fgMinuitStepSize = 0.1;                 // (usually not needed to be changed) step size in minimization
50 Bool_t AliUnfolding::fgMinimumInitialValue = kFALSE;          // set all initial values at least to the smallest value among the initial values
51 Float_t AliUnfolding::fgMinimumInitialValueFix = -1;
52 Bool_t AliUnfolding::fgNormalizeInput = kFALSE;                  // normalize input spectrum
53 Float_t AliUnfolding::fgNotFoundEvents = 0;
54 Bool_t AliUnfolding::fgSkipBin0InChi2 = kFALSE;
55
56 Float_t AliUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
57 Int_t   AliUnfolding::fgBayesianIterations = 10;          // number of iterations in Bayesian method
58
59 Bool_t AliUnfolding::fgDebug = kFALSE;
60
61 Int_t AliUnfolding::fgCallCount = 0;
62
63 ClassImp(AliUnfolding)
64
65 //____________________________________________________________________
66 void AliUnfolding::SetUnfoldingMethod(MethodType methodType)
67 {
68   // set unfolding method
69   fgMethodType = methodType; 
70   
71   const char* name = 0;
72   switch (methodType)
73   {
74     case kInvalid: name = "INVALID"; break;
75     case kChi2Minimization: name = "Chi2 Minimization"; break;
76     case kBayesian: name = "Bayesian unfolding"; break;
77     case kFunction: name = "Functional fit"; break;
78   }
79   Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name);
80 }
81
82 //____________________________________________________________________
83 void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit) 
84
85   // enable the creation of a overflow bin that includes all statistics below the given limit
86   
87   fgOverflowBinLimit = overflowBinLimit; 
88   
89   Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit);
90 }
91
92 //____________________________________________________________________
93 void AliUnfolding::SetSkipBinsBegin(Int_t nBins)
94 {
95   // set number of skipped bins in regularization
96   
97   fgSkipBinsBegin = nBins;
98   
99   Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin);
100 }
101
102 //____________________________________________________________________
103 void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded) 
104
105   // set number of bins in the input (measured) distribution and in the unfolded distribution
106   fgMaxInput = nMeasured; 
107   fgMaxParams = nUnfolded; 
108   
109   if (fgCorrelationMatrix)
110   {
111     delete fgCorrelationMatrix;
112     fgCorrelationMatrix = 0;
113   }
114   if (fgCorrelationMatrixSquared)
115   {
116     fgCorrelationMatrixSquared = 0;
117     delete fgCorrelationMatrixSquared;
118   }
119   if (fgCorrelationCovarianceMatrix)
120   {
121     delete fgCorrelationCovarianceMatrix;
122     fgCorrelationCovarianceMatrix = 0;
123   }
124   if (fgCurrentESDVector)
125   {
126     delete fgCurrentESDVector;
127     fgCurrentESDVector = 0;
128   }
129   if (fgEntropyAPriori)
130   {
131     delete fgEntropyAPriori;
132     fgEntropyAPriori = 0;
133   }
134   if (fgEfficiency)
135   {
136     delete fgEfficiency;
137     fgEfficiency = 0;
138   }
139   if (fgBinWidths)
140   {
141     delete fgBinWidths;
142     fgBinWidths = 0;
143   }
144   
145   Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
146 }
147
148 //____________________________________________________________________
149 void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight)
150 {
151   //
152   // sets the parameters for chi2 minimization
153   //
154
155   fgRegularizationType = type;
156   fgRegularizationWeight = weight;
157
158   Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight);
159 }
160
161 //____________________________________________________________________
162 void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
163 {
164   //
165   // sets the parameters for Bayesian unfolding
166   //
167
168   fgBayesianSmoothing = smoothing;
169   fgBayesianIterations = nIterations;
170
171   Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
172 }
173
174 //____________________________________________________________________
175 void AliUnfolding::SetFunction(TF1* function)
176 {
177   // set function for unfolding with a fit function
178   
179   fgFitFunction = function;
180   
181   Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar());
182 }
183
184 //____________________________________________________________________
185 Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
186 {
187   // unfolds with unfolding method fgMethodType
188   //
189   // parameters:
190   //  correlation: response matrix as measured vs. generated
191   //  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.
192   //  measured:    the measured spectrum
193   //  initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions.
194   //  result:      target for the unfolded result
195   //  check:       depends on the unfolding method, see comments in specific functions
196   //
197   //  return code: see UnfoldWithMinuit/UnfoldWithBayesian/UnfoldWithFunction
198
199   if (fgMaxInput == -1)
200   {
201     Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
202     fgMaxInput = measured->GetNbinsX();
203   }
204   if (fgMaxParams == -1)
205   {
206     Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
207     fgMaxParams = measured->GetNbinsX();
208   }
209
210   if (fgOverflowBinLimit > 0)
211     CreateOverflowBin(correlation, measured);
212     
213   switch (fgMethodType)
214   {
215     case kInvalid:
216     {
217       Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting...");
218       return -1;
219     }
220     case kChi2Minimization:
221       return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check);
222     case kBayesian:
223       return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result);
224     case kFunction:
225       return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result);
226   }
227   return -1;
228 }
229
230 //____________________________________________________________________
231 void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency)
232 {
233   // fill static variables needed for minuit fit
234
235   if (!fgCorrelationMatrix)
236     fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
237   if (!fgCorrelationMatrixSquared)
238     fgCorrelationMatrixSquared = new TMatrixD(fgMaxInput, fgMaxParams);
239   if (!fgCorrelationCovarianceMatrix)
240     fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
241   if (!fgCurrentESDVector)
242     fgCurrentESDVector = new TVectorD(fgMaxInput);
243   if (!fgEntropyAPriori)
244     fgEntropyAPriori = new TVectorD(fgMaxParams);
245   if (!fgEfficiency)
246     fgEfficiency = new TVectorD(fgMaxParams);
247   if (!fgBinWidths)
248     fgBinWidths = new TVectorD(fgMaxParams);
249     
250   fgCorrelationMatrix->Zero();
251   fgCorrelationCovarianceMatrix->Zero();
252   fgCurrentESDVector->Zero();
253   fgEntropyAPriori->Zero();
254
255   // normalize correction for given nPart
256   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
257   {
258     Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
259     if (sum <= 0)
260       continue;
261     Float_t maxValue = 0;
262     Int_t maxBin = -1;
263     for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
264     {
265       // find most probably value
266       if (maxValue < correlation->GetBinContent(i, j))
267       {
268         maxValue = correlation->GetBinContent(i, j);
269         maxBin = j;
270       }
271
272       // npart sum to 1
273       correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i));
274       correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
275
276       if (i <= fgMaxParams && j <= fgMaxInput)
277       {
278         (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
279         (*fgCorrelationMatrixSquared)(j-1, i-1) = correlation->GetBinContent(i, j) * correlation->GetBinContent(i, j);
280       }
281     }
282
283     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
284   }
285     
286   //normalize measured
287   Float_t smallestError = 1;
288   if (fgNormalizeInput)
289   {
290     Float_t sumMeasured = measured->Integral();
291     measured->Scale(1.0 / sumMeasured);
292     smallestError /= sumMeasured;
293   }
294   
295   for (Int_t i=0; i<fgMaxInput; ++i)
296   {
297     (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
298     if (measured->GetBinError(i+1) > 0)
299     {
300       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
301     }
302     else // in this case put error of 1, otherwise 0 bins are not added to the chi2...
303       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / smallestError / smallestError;
304
305     if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
306       (*fgCorrelationCovarianceMatrix)(i, i) = 0;
307     //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
308   }
309
310   // efficiency is expected to match bin width of result
311   for (Int_t i=0; i<fgMaxParams; ++i)
312   {
313     (*fgEfficiency)(i) = efficiency->GetBinContent(i+1);
314     (*fgBinWidths)(i) = efficiency->GetXaxis()->GetBinWidth(i+1);
315   }
316 }
317
318 //____________________________________________________________________
319 Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
320 {
321   //
322   // implementation of unfolding (internal function)
323   //
324   // unfolds <measured> using response from <correlation> and effiency <efficiency>
325   // output is in <result>
326   // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
327   //   negative values in initialConditions mean that the given parameter is fixed to the absolute of the value
328   // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
329   //
330   // returns minuit status (0 = success), (-1 when check was set)
331   //
332
333   SetStaticVariables(correlation, measured, efficiency);
334   
335   // Initialize TMinuit via generic fitter interface
336   Int_t params = fgMaxParams;
337   if (fgNotFoundEvents > 0)
338     params++;
339     
340   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, params);
341   Double_t arglist[100];
342
343   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
344   arglist[0] = 0;
345   minuit->ExecuteCommand("SET PRINT", arglist, 1);
346
347   // however, enable warnings
348   //minuit->ExecuteCommand("SET WAR", arglist, 0);
349
350   // set minimization function
351   minuit->SetFCN(Chi2Function);
352
353   for (Int_t i=0; i<fgMaxParams; i++)
354     (*fgEntropyAPriori)[i] = 1;
355
356   // set initial conditions as a-priori distribution for MRX regularization
357   /*
358   for (Int_t i=0; i<fgMaxParams; i++)
359     if (initialConditions && initialConditions->GetBinContent(i+1) > 0)
360       (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
361   */
362
363   if (!initialConditions) {
364     initialConditions = measured;
365   } else {
366     Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions...");
367     //new TCanvas; initialConditions->DrawCopy();
368     if (fgNormalizeInput)
369       initialConditions->Scale(1.0 / initialConditions->Integral());
370   }
371
372   // extract minimum value from initial conditions (if we set a value to 0 it will stay 0)
373   Float_t minValue = 1e100;
374   if (fgMinimumInitialValueFix < 0)
375   {
376     for (Int_t i=0; i<fgMaxParams; ++i)
377     {
378       Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
379       if (initialConditions->GetBinContent(bin) > 0)
380         minValue = TMath::Min(minValue, (Float_t) initialConditions->GetBinContent(bin));
381     }
382   }
383   else
384     minValue = fgMinimumInitialValueFix;
385   
386   Double_t* results = new Double_t[fgMaxParams+1];
387   for (Int_t i=0; i<fgMaxParams; ++i)
388   {
389     Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
390     results[i] = initialConditions->GetBinContent(bin);
391
392     Bool_t fix = kFALSE;
393     if (results[i] < 0)
394     {
395       fix = kTRUE;
396       results[i] = -results[i];
397     }
398  
399     if (!fix && fgMinimumInitialValue && results[i] < minValue)
400       results[i] = minValue;
401       
402     // minuit sees squared values to prevent it from going negative...
403     results[i] = TMath::Sqrt(results[i]);
404
405     minuit->SetParameter(i, Form("param%d", i), results[i], (fix) ? 0 : fgMinuitStepSize, 0, 0);
406   }
407   if (fgNotFoundEvents > 0)
408   {
409     results[fgMaxParams] = efficiency->GetBinContent(1);
410     minuit->SetParameter(fgMaxParams, "vtx0", results[fgMaxParams], fgMinuitStepSize / 100, 0.01, 0.80);
411   }
412   
413   Int_t dummy = 0;
414   Double_t chi2 = 0;
415   Chi2Function(dummy, 0, chi2, results, 0);
416   printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2);
417
418   if (check)
419   {
420     DrawGuess(results);
421     delete[] results;
422     return -1;
423   }
424
425   // first param is number of iterations, second is precision....
426   arglist[0] = 1e6;
427   //arglist[1] = 1e-5;
428   //minuit->ExecuteCommand("SCAN", arglist, 0);
429   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
430   Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status);
431   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
432   //minuit->ExecuteCommand("MINOS", arglist, 0);
433
434   if (fgNotFoundEvents > 0)
435   {
436     results[fgMaxParams] = minuit->GetParameter(fgMaxParams);
437     Printf("Efficiency for bin 0 changed from %f to %f", efficiency->GetBinContent(1), results[fgMaxParams]);
438     efficiency->SetBinContent(1, results[fgMaxParams]);
439   }
440   
441   for (Int_t i=0; i<fgMaxParams; ++i)
442   {
443     results[i] = minuit->GetParameter(i);
444     Double_t value = results[i] * results[i];
445     // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
446     Double_t error = minuit->GetParError(i) * results[i];
447     
448     if (efficiency)
449     {   
450       if (efficiency->GetBinContent(i+1) > 0)
451       {
452         value /= efficiency->GetBinContent(i+1);
453         error /= efficiency->GetBinContent(i+1);
454       }
455       else
456       {
457         value = 0;
458         error = 0;
459       }
460     }
461     
462     result->SetBinContent(i+1, value);
463     result->SetBinError(i+1, error);
464   }
465   
466   fgCallCount = 0;
467   Chi2Function(dummy, 0, chi2, results, 0);
468   printf("AliUnfolding::UnfoldWithMinuit: Chi2 of final parameters is = %f\n", chi2);
469   
470   if (fgDebug)
471     DrawGuess(results);
472
473   delete[] results;
474
475   return status;
476 }
477
478 //____________________________________________________________________
479 Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult)
480 {
481   //
482   // unfolds a spectrum using the Bayesian method
483   //
484   
485   if (measured->Integral() <= 0)
486   {
487     Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
488     return -1;
489   }
490
491   const Int_t kStartBin = 0;
492
493   Int_t kMaxM = fgMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
494   Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
495
496   // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
497   const Double_t kConvergenceLimit = kMaxT * 1e-6;
498
499   // store information in arrays, to increase processing speed (~ factor 5)
500   Double_t* measuredCopy = new Double_t[kMaxM];
501   Double_t* measuredError = new Double_t[kMaxM];
502   Double_t* prior = new Double_t[kMaxT];
503   Double_t* result = new Double_t[kMaxT];
504   Double_t* efficiency = new Double_t[kMaxT];
505   Double_t* binWidths = new Double_t[kMaxT];
506
507   Double_t** response = new Double_t*[kMaxT];
508   Double_t** inverseResponse = new Double_t*[kMaxT];
509   for (Int_t i=0; i<kMaxT; i++)
510   {
511     response[i] = new Double_t[kMaxM];
512     inverseResponse[i] = new Double_t[kMaxM];
513   }
514
515   // for normalization
516   Float_t measuredIntegral = measured->Integral();
517   for (Int_t m=0; m<kMaxM; m++)
518   {
519     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
520     measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
521
522     for (Int_t t=0; t<kMaxT; t++)
523     {
524       response[t][m] = correlation->GetBinContent(t+1, m+1);
525       inverseResponse[t][m] = 0;
526     }
527   }
528
529   for (Int_t t=0; t<kMaxT; t++)
530   {
531     if (aEfficiency)
532     {
533       efficiency[t] = aEfficiency->GetBinContent(t+1);
534     }
535     else
536       efficiency[t] = 1;
537       
538     prior[t] = measuredCopy[t];
539     result[t] = 0;
540     binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1);
541   }
542
543   // pick prior distribution
544   if (initialConditions)
545   {
546     printf("Using different starting conditions...\n");
547     // for normalization
548     Float_t inputDistIntegral = initialConditions->Integral();
549     for (Int_t i=0; i<kMaxT; i++)
550       prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
551   }
552
553   //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
554   
555   //new TCanvas;
556   // unfold...
557   for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations < 0; i++)
558   {
559     if (fgDebug)
560       Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
561
562     // calculate IR from Bayes theorem
563     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
564
565     Double_t chi2Measured = 0;
566     for (Int_t m=0; m<kMaxM; m++)
567     {
568       Float_t norm = 0;
569       for (Int_t t = kStartBin; t<kMaxT; t++)
570         norm += response[t][m] * prior[t];
571
572       // calc. chi2: (measured - response * prior) / error
573       if (measuredError[m] > 0)
574       {
575         Double_t value = (measuredCopy[m] - norm) / measuredError[m];
576         chi2Measured += value * value;
577       }
578
579       if (norm > 0)
580       {
581         for (Int_t t = kStartBin; t<kMaxT; t++)
582           inverseResponse[t][m] = response[t][m] * prior[t] / norm;
583       }
584       else
585       {
586         for (Int_t t = kStartBin; t<kMaxT; t++)
587           inverseResponse[t][m] = 0;
588       }
589     }
590     //Printf("chi2Measured of the last prior is %e", chi2Measured);
591
592     for (Int_t t = kStartBin; t<kMaxT; t++)
593     {
594       Float_t value = 0;
595       for (Int_t m=0; m<kMaxM; m++)
596         value += inverseResponse[t][m] * measuredCopy[m];
597
598       if (efficiency[t] > 0)
599         result[t] = value / efficiency[t];
600       else
601         result[t] = 0;
602     }
603     
604     /* 
605     // draw intermediate result
606     for (Int_t t=0; t<kMaxT; t++)
607     {
608       aResult->SetBinContent(t+1, result[t]);
609     }
610     aResult->SetMarkerStyle(24+i);
611     aResult->SetMarkerColor(2);
612     aResult->DrawCopy((i == 0) ? "P" : "PSAME");
613     */
614  
615     Double_t chi2LastIter = 0;
616     // regularization (simple smoothing)
617     for (Int_t t=kStartBin; t<kMaxT; t++)
618     {
619       Float_t newValue = 0;
620       
621       // 0 bin excluded from smoothing
622       if (t > kStartBin+2 && t<kMaxT-1)
623       {
624         Float_t average = (result[t-1] / binWidths[t-1] + result[t] / binWidths[t] + result[t+1] / binWidths[t+1]) / 3 * binWidths[t];
625
626         // weight the average with the regularization parameter
627         newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
628       }
629       else
630         newValue = result[t];
631
632       // calculate chi2 (change from last iteration)
633       if (prior[t] > 1e-5)
634       {
635         Double_t diff = (prior[t] - newValue) / prior[t];
636         chi2LastIter += diff * diff;
637       }
638
639       prior[t] = newValue;
640     }
641     //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
642     //convergence->Fill(i+1, chi2LastIter);
643
644     if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit)
645     {
646       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);
647       break;
648     }
649   } // end of iterations
650
651   //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
652   //delete convergence;
653
654   Float_t factor = 1;
655   if (!fgNormalizeInput)
656     factor = measuredIntegral;
657   for (Int_t t=0; t<kMaxT; t++)
658     aResult->SetBinContent(t+1, result[t] * factor);
659
660   delete[] measuredCopy;
661   delete[] measuredError;
662   delete[] prior;
663   delete[] result;
664   delete[] efficiency;
665   delete[] binWidths;
666
667   for (Int_t i=0; i<kMaxT; i++)
668   {
669     delete[] response[i];
670     delete[] inverseResponse[i];
671   }
672   delete[] response;
673   delete[] inverseResponse;
674   
675   return 0;
676
677   // ********
678   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
679
680   /*printf("Calculating covariance matrix. This may take some time...\n");
681
682   // check if this is the right one...
683   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
684
685   Int_t xBins = hInverseResponseBayes->GetNbinsX();
686   Int_t yBins = hInverseResponseBayes->GetNbinsY();
687
688   // calculate "unfolding matrix" Mij
689   Float_t matrixM[251][251];
690   for (Int_t i=1; i<=xBins; i++)
691   {
692     for (Int_t j=1; j<=yBins; j++)
693     {
694       if (fCurrentEfficiency->GetBinContent(i) > 0)
695         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
696       else
697         matrixM[i-1][j-1] = 0;
698     }
699   }
700
701   Float_t* vectorn = new Float_t[yBins];
702   for (Int_t j=1; j<=yBins; j++)
703     vectorn[j-1] = fCurrentESD->GetBinContent(j);
704
705   // first part of covariance matrix, depends on input distribution n(E)
706   Float_t cov1[251][251];
707
708   Float_t nEvents = fCurrentESD->Integral(); // N
709
710   xBins = 20;
711   yBins = 20;
712
713   for (Int_t k=0; k<xBins; k++)
714   {
715     printf("In Cov1: %d\n", k);
716     for (Int_t l=0; l<yBins; l++)
717     {
718       cov1[k][l] = 0;
719
720       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
721       for (Int_t j=0; j<yBins; j++)
722         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
723           * (1.0 - vectorn[j] / nEvents);
724
725       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
726       for (Int_t i=0; i<yBins; i++)
727         for (Int_t j=0; j<yBins; j++)
728         {
729           if (i == j)
730             continue;
731           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
732             * vectorn[j] / nEvents;
733          }
734     }
735   }
736
737   printf("Cov1 finished\n");
738
739   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
740   cov->Reset();
741
742   for (Int_t i=1; i<=xBins; i++)
743     for (Int_t j=1; j<=yBins; j++)
744       cov->SetBinContent(i, j, cov1[i-1][j-1]);
745
746   new TCanvas;
747   cov->Draw("COLZ");
748
749   // second part of covariance matrix, depends on response matrix
750   Float_t cov2[251][251];
751
752   // Cov[P(Er|Cu), P(Es|Cu)] term
753   Float_t covTerm[100][100][100];
754   for (Int_t r=0; r<yBins; r++)
755     for (Int_t u=0; u<xBins; u++)
756       for (Int_t s=0; s<yBins; s++)
757       {
758         if (r == s)
759           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
760             * (1.0 - hResponse->GetBinContent(u+1, r+1));
761         else
762           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
763             * hResponse->GetBinContent(u+1, s+1);
764       }
765
766   for (Int_t k=0; k<xBins; k++)
767     for (Int_t l=0; l<yBins; l++)
768     {
769       cov2[k][l] = 0;
770       printf("In Cov2: %d %d\n", k, l);
771       for (Int_t i=0; i<yBins; i++)
772         for (Int_t j=0; j<yBins; j++)
773         {
774           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
775           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
776           Float_t tmpCov = 0;
777           for (Int_t r=0; r<yBins; r++)
778             for (Int_t u=0; u<xBins; u++)
779               for (Int_t s=0; s<yBins; s++)
780               {
781                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
782                   || hResponse->GetBinContent(u+1, i+1) == 0)
783                   continue;
784
785                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
786                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
787                   * covTerm[r][u][s];
788               }
789
790           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
791         }
792     }
793
794   printf("Cov2 finished\n");
795
796   for (Int_t i=1; i<=xBins; i++)
797     for (Int_t j=1; j<=yBins; j++)
798       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
799
800   new TCanvas;
801   cov->Draw("COLZ");*/
802 }
803
804 //____________________________________________________________________
805 Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
806 {
807   // homogenity term for minuit fitting
808   // pure function of the parameters
809   // prefers constant function (pol0)
810
811   Double_t chi2 = 0;
812
813   for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
814   {
815     Double_t right  = params[i] / (*fgBinWidths)(i);
816     Double_t left   = params[i-1] / (*fgBinWidths)(i-1);
817
818     if (left != 0)
819     {
820       Double_t diff = (right - left);
821       chi2 += diff * diff / left / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
822     }
823   }
824
825   return chi2 / 100.0;
826 }
827
828 //____________________________________________________________________
829 Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
830 {
831   // homogenity term for minuit fitting
832   // pure function of the parameters
833   // prefers linear function (pol1)
834
835   Double_t chi2 = 0;
836
837   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
838   {
839     if (params[i-1] == 0)
840       continue;
841
842     Double_t right  = params[i] / (*fgBinWidths)(i);
843     Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
844     Double_t left   = params[i-2] / (*fgBinWidths)(i-2);
845
846     Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
847     Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2);
848
849     //Double_t diff = (der1 - der2) / middle;
850     //chi2 += diff * diff;
851     chi2 += (der1 - der2) * (der1 - der2) / middle * (*fgBinWidths)(i-1);
852   }
853
854   return chi2;
855 }
856
857 //____________________________________________________________________
858 Double_t AliUnfolding::RegularizationLog(TVectorD& params)
859 {
860   // homogenity term for minuit fitting
861   // pure function of the parameters
862   // prefers linear function (pol1)
863
864   Double_t chi2 = 0;
865
866   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
867   {
868     if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
869      continue;
870
871     Double_t right  = log(params[i] / (*fgBinWidths)(i));
872     Double_t middle = log(params[i-1] / (*fgBinWidths)(i-1));
873     Double_t left   = log(params[i-2] / (*fgBinWidths)(i-2));
874     
875     Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
876     Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2);
877     
878     //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
879
880     //if (fgCallCount == 0)
881     //  Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
882     chi2 += (der1 - der2) * (der1 - der2);// / error;
883   }
884
885   return chi2;
886 }
887
888 //____________________________________________________________________
889 Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
890 {
891   // homogenity term for minuit fitting
892   // pure function of the parameters
893   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
894   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
895
896   Double_t chi2 = 0;
897
898   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
899   {
900     Double_t right  = params[i];
901     Double_t middle = params[i-1];
902     Double_t left   = params[i-2];
903
904     Double_t der1 = (right - middle);
905     Double_t der2 = (middle - left);
906
907     Double_t diff = (der1 - der2);
908
909     chi2 += diff * diff;
910   }
911
912   return chi2 * 1e4;
913 }
914
915 //____________________________________________________________________
916 Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
917 {
918   // homogenity term for minuit fitting
919   // pure function of the parameters
920   // calculates entropy, from
921   // The method of reduced cross-entropy (M. Schmelling 1993)
922
923   Double_t paramSum = 0;
924   
925   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
926     paramSum += params[i];
927
928   Double_t chi2 = 0;
929   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
930   {
931     Double_t tmp = params[i] / paramSum;
932     //Double_t tmp = params[i];
933     if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
934     {
935       chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
936     }
937     else
938       chi2 += 100;
939   }
940
941   return -chi2;
942 }
943
944 //____________________________________________________________________
945 Double_t AliUnfolding::RegularizationRatio(TVectorD& params)
946 {
947   // homogenity term for minuit fitting
948   // pure function of the parameters
949
950   Double_t chi2 = 0;
951
952   for (Int_t i=5+fgSkipBinsBegin; i<fgMaxParams; ++i)
953   {
954     if (params[i-1] == 0 || params[i] == 0)
955       continue;
956
957     Double_t right  = params[i] / (*fgBinWidths)(i);
958     Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
959     Double_t left   = params[i-2] / (*fgBinWidths)(i-2);
960     Double_t left2   = params[i-3] / (*fgBinWidths)(i-3);
961     Double_t left3   = params[i-4] / (*fgBinWidths)(i-4);
962     Double_t left4   = params[i-5] / (*fgBinWidths)(i-5);
963
964     //Double_t diff = left / middle - middle / right;
965     //Double_t diff = 2 * left / middle - middle / right - left2 / left;
966     Double_t diff = 4 * left2 / left - middle / right - left / middle - left3 / left2 - left4 / left3;
967     
968     chi2 += diff * diff;// / middle;
969   }
970
971   return chi2;
972 }
973
974 //____________________________________________________________________
975 void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
976 {
977   //
978   // fit function for minuit
979   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
980   //
981   
982   // TODO use static members for the variables here to speed up processing (no construction/deconstruction)
983
984   // d
985   TVectorD paramsVector(fgMaxParams);
986   for (Int_t i=0; i<fgMaxParams; ++i)
987     paramsVector[i] = params[i] * params[i];
988
989   // calculate penalty factor
990   Double_t penaltyVal = 0;
991   switch (fgRegularizationType)
992   {
993     case kNone:       break;
994     case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
995     case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
996     case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
997     case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
998     case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
999     case kRatio:      penaltyVal = RegularizationRatio(paramsVector); break;
1000   }
1001
1002   //if (penaltyVal > 1e10)
1003   //  paramsVector2.Print();
1004
1005   penaltyVal *= fgRegularizationWeight;
1006
1007   // Ad
1008   TVectorD measGuessVector(fgMaxInput);
1009   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1010
1011   // Ad - m
1012   measGuessVector -= (*fgCurrentESDVector);
1013   
1014 #if 0
1015   // new error calcuation using error on the guess
1016   
1017   // error from guess
1018   TVectorD errorGuessVector(fgMaxInput);
1019   //errorGuessVector = (*fgCorrelationMatrixSquared) * paramsVector;
1020   errorGuessVector = (*fgCorrelationMatrix) * paramsVector;
1021
1022   Double_t chi2FromFit = 0;
1023   for (Int_t i=0; i<fgMaxInput; ++i)
1024   {
1025     Float_t error = 1;
1026     if (errorGuessVector(i) > 0)
1027       error = errorGuessVector(i);
1028     chi2FromFit += measGuessVector(i) * measGuessVector(i) / error;
1029   }
1030   //measGuessVector.Print();
1031
1032 #else
1033   // old error calcuation using the error on the measured
1034   TVectorD copy(measGuessVector);
1035
1036   // (Ad - m) W
1037   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1038   // normal way is like this:
1039   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1040   // optimized way like this:
1041   for (Int_t i=0; i<fgMaxInput; ++i)
1042     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1043
1044   //measGuessVector.Print();
1045
1046   if (fgSkipBin0InChi2)
1047     measGuessVector[0] = 0;
1048
1049   // (Ad - m) W (Ad - m)
1050   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1051   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see UnfoldWithMinuit)
1052   Double_t chi2FromFit = measGuessVector * copy * 1e6;
1053 #endif
1054
1055   Double_t notFoundEventsConstraint = 0;
1056   Double_t currentNotFoundEvents = 0;
1057   Double_t errorNotFoundEvents = 0;
1058   
1059   if (fgNotFoundEvents > 0)
1060   {
1061     for (Int_t i=0; i<fgMaxParams; ++i)
1062     {
1063       Double_t eff = (1.0 / (*fgEfficiency)(i) - 1);
1064       if (i == 0)
1065         eff = (1.0 / params[fgMaxParams] - 1);
1066       currentNotFoundEvents += eff * paramsVector(i);
1067       errorNotFoundEvents += eff * eff * paramsVector(i); // error due to guess (paramsVector)
1068       if (i <= 3)
1069         errorNotFoundEvents += (eff * 0.03) * (eff * 0.03) * paramsVector(i) * paramsVector(i); // error on eff
1070       //if ((fgCallCount % 10000) == 0)
1071       //  Printf("%d %f %f %f", i, (*fgEfficiency)(i), paramsVector(i), currentNotFoundEvents);
1072     }
1073     errorNotFoundEvents += fgNotFoundEvents;
1074     // TODO add error on background, background estimate
1075     
1076     notFoundEventsConstraint = (currentNotFoundEvents - fgNotFoundEvents) * (currentNotFoundEvents - fgNotFoundEvents) / errorNotFoundEvents;
1077   }
1078   
1079   Float_t avg = 0;
1080   Float_t sum = 0;
1081   Float_t currentMult = 0;
1082   // try to match dn/deta
1083   for (Int_t i=0; i<fgMaxParams; ++i)
1084   {
1085     avg += paramsVector(i) * currentMult;
1086     sum += paramsVector(i);
1087     currentMult += (*fgBinWidths)(i);
1088   }
1089   avg /= sum;
1090   Float_t chi2Avg = 0; //(avg - 3.73) * (avg - 3.73) * 100;
1091
1092   chi2 = chi2FromFit + penaltyVal + notFoundEventsConstraint + chi2Avg;
1093   
1094   if ((fgCallCount++ % 1000) == 0)
1095   {
1096     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], avg, chi2FromFit, penaltyVal, notFoundEventsConstraint, chi2Avg, chi2);
1097     //for (Int_t i=0; i<fgMaxInput; ++i)
1098     //  Printf("%d: %f", i, measGuessVector(i) * copy(i) * 1e6);
1099   }
1100 }
1101
1102 //____________________________________________________________________
1103 void AliUnfolding::DrawGuess(Double_t *params)
1104 {
1105   //
1106   // draws residuals of solution suggested by params and effect of regularization
1107   //
1108
1109   // d
1110   TVectorD paramsVector(fgMaxParams);
1111   for (Int_t i=0; i<fgMaxParams; ++i)
1112     paramsVector[i] = params[i] * params[i];
1113
1114   // Ad
1115   TVectorD measGuessVector(fgMaxInput);
1116   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1117
1118   // Ad - m
1119   measGuessVector -= (*fgCurrentESDVector);
1120
1121   TVectorD copy(measGuessVector);
1122   //copy.Print();
1123
1124   // (Ad - m) W
1125   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1126   // normal way is like this:
1127   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1128   // optimized way like this:
1129   for (Int_t i=0; i<fgMaxInput; ++i)
1130     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1131
1132   // (Ad - m) W (Ad - m)
1133   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1134   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
1135   //Double_t chi2FromFit = measGuessVector * copy * 1e6;
1136
1137   // draw residuals
1138   TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
1139   for (Int_t i=0; i<fgMaxInput; ++i)
1140     residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
1141   new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
1142
1143   // draw penalty
1144   TH1* penalty = GetPenaltyPlot(params);
1145   
1146   new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
1147 }
1148
1149 //____________________________________________________________________
1150 TH1* AliUnfolding::GetPenaltyPlot(TH1* corrected)
1151 {
1152   // draws the penalty factors as function of multiplicity of the current selected regularization
1153
1154   Double_t* params = new Double_t[fgMaxParams];
1155   for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
1156     params[i] = corrected->GetBinContent(i+1);
1157   
1158   TH1* penalty = GetPenaltyPlot(params);
1159   
1160   delete[] params;
1161   
1162   return penalty;
1163 }
1164
1165 //____________________________________________________________________
1166 TH1* AliUnfolding::GetPenaltyPlot(Double_t* params)
1167 {
1168   // draws the penalty factors as function of multiplicity of the current selected regularization
1169   
1170   TH1* penalty = new TH1F("penalty", ";unfolded multiplicity;penalty factor", fgMaxParams, -0.5, fgMaxParams - 0.5);
1171
1172   for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
1173   {
1174     Double_t diff = 0;
1175     if (fgRegularizationType == kPol0)
1176     {
1177       Double_t right  = params[i] / (*fgBinWidths)(i);
1178       Double_t left   = params[i-1] / (*fgBinWidths)(i-1);
1179
1180       if (left != 0)
1181       {
1182         Double_t diffTmp = (right - left);
1183         diff = diffTmp * diffTmp / left / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2) / 100;
1184       }
1185     } 
1186     if (fgRegularizationType == kPol1 && i > 1+fgSkipBinsBegin) 
1187     {
1188       if (params[i-1] == 0)
1189         continue;
1190
1191       Double_t right  = params[i] / (*fgBinWidths)(i);
1192       Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
1193       Double_t left   = params[i-2] / (*fgBinWidths)(i-2);
1194
1195       Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
1196       Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2);
1197
1198       diff = (der1 - der2) * (der1 - der2) / middle;
1199     }
1200     if (fgRegularizationType == kLog && i > 1+fgSkipBinsBegin) 
1201     {
1202       if (params[i-1] == 0)
1203         continue;
1204   
1205       Double_t right  = log(params[i]);
1206       Double_t middle = log(params[i-1]);
1207       Double_t left   = log(params[i-2]);
1208       
1209       Double_t der1 = (right - middle);
1210       Double_t der2 = (middle - left);
1211   
1212       //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
1213       //Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
1214        
1215       diff = (der1 - der2) * (der1 - der2);// / error;
1216     }
1217
1218     penalty->SetBinContent(i, diff);
1219     
1220     //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
1221   }
1222   
1223   return penalty;
1224 }
1225
1226 //____________________________________________________________________
1227 void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
1228 {
1229   //
1230   // fit function for minuit
1231   // uses the TF1 stored in fgFitFunction
1232   //
1233
1234   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1235     fgFitFunction->SetParameter(i, params[i]);
1236
1237   Double_t* params2 = new Double_t[fgMaxParams];
1238
1239   for (Int_t i=0; i<fgMaxParams; ++i)
1240     params2[i] = fgFitFunction->Eval(i);
1241
1242   Chi2Function(unused1, unused2, chi2, params2, unused3);
1243   
1244   delete[] params2;
1245
1246   if (fgDebug)
1247     Printf("%f", chi2);
1248 }
1249
1250 //____________________________________________________________________
1251 Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult)
1252 {
1253   //
1254   // correct spectrum using minuit chi2 method applying a functional fit
1255   //
1256   
1257   if (!fgFitFunction)
1258   {
1259     Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting.");
1260     return -1;
1261   }    
1262   
1263   SetChi2Regularization(kNone, 0);
1264   
1265   SetStaticVariables(correlation, measured, efficiency);
1266   
1267   // Initialize TMinuit via generic fitter interface
1268   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
1269
1270   minuit->SetFCN(TF1Function);
1271   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1272   {
1273     Double_t lower, upper;
1274     fgFitFunction->GetParLimits(i, lower, upper);
1275     minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper);
1276   }
1277
1278   Double_t arglist[100];
1279   arglist[0] = 0;
1280   minuit->ExecuteCommand("SET PRINT", arglist, 1);
1281   minuit->ExecuteCommand("MIGRAD", arglist, 0);
1282   //minuit->ExecuteCommand("MINOS", arglist, 0);
1283
1284   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1285     fgFitFunction->SetParameter(i, minuit->GetParameter(i));
1286
1287   for (Int_t i=0; i<fgMaxParams; ++i)
1288   {
1289     Double_t value = fgFitFunction->Eval(i);
1290     if (fgDebug)
1291       Printf("%d : %f", i, value);
1292     
1293     if (efficiency)
1294     {
1295       if (efficiency->GetBinContent(i+1) > 0)
1296       {
1297         value /= efficiency->GetBinContent(i+1);
1298       }
1299       else
1300         value = 0;
1301     }
1302     aResult->SetBinContent(i+1, value);
1303     aResult->SetBinError(i+1, 0);
1304   }
1305   
1306   return 0;
1307 }
1308
1309 //____________________________________________________________________
1310 void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
1311 {
1312   // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins
1313   // The same limit is applied to the correlation  
1314   
1315   Int_t lastBin = 0;
1316   for (Int_t i=1; i<=measured->GetNbinsX(); ++i)
1317   {
1318     if (measured->GetBinContent(i) <= fgOverflowBinLimit)
1319     {
1320       lastBin = i;
1321       break;
1322     }
1323   }
1324   
1325   if (lastBin == 0)
1326   {
1327     Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit);
1328     return;
1329   }
1330   
1331   Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin);
1332   
1333   TCanvas* canvas = 0;
1334
1335   if (fgDebug)
1336   {
1337     canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
1338     canvas->Divide(2, 2);
1339
1340     canvas->cd(1);
1341     measured->SetStats(kFALSE);
1342     measured->DrawCopy();
1343     gPad->SetLogy();
1344
1345     canvas->cd(2);
1346     correlation->SetStats(kFALSE);
1347     correlation->DrawCopy("COLZ");
1348   }
1349
1350   measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX()));
1351   for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i)
1352   {
1353     measured->SetBinContent(i, 0);
1354     measured->SetBinError(i, 0);
1355   }
1356   // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1357   measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin)));
1358
1359   Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin));
1360
1361   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1362   {
1363     correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY()));
1364     // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1365     correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin)));
1366
1367     for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j)
1368     {
1369       correlation->SetBinContent(i, j, 0);
1370       correlation->SetBinError(i, j, 0);
1371     }
1372   }
1373
1374   Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!");
1375
1376   if (canvas)
1377   {
1378     canvas->cd(3);
1379     measured->DrawCopy();
1380     gPad->SetLogy();
1381
1382     canvas->cd(4);
1383     correlation->DrawCopy("COLZ");
1384   }
1385 }
1386
1387 Int_t AliUnfolding::UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result)
1388 {
1389   // unfolds and assigns bias as errors with Eq. 19 of Cowan, "a survey of unfolding methods for particle physics"
1390   // b_i = sum_j dmu_i/dn_j (nu_j - n_j) with nu_j as folded guess, n_j as data
1391   // dmu_i/dn_j is found numerically by changing the bin content and re-unfolding
1392   //
1393   // return code: 0 (success) -1 (error: from Unfold(...) )
1394
1395   if (Unfold(correlation, efficiency, measured, initialConditions, result) != 0)
1396     return -1;
1397
1398   TMatrixD matrix(fgMaxInput, fgMaxParams);
1399   
1400   TH1* newResult[4];
1401   for (Int_t loop=0; loop<4; loop++)
1402     newResult[loop] = (TH1*) result->Clone(Form("newresult_%d", loop));
1403     
1404   // change bin-by-bin and built matrix of effects
1405   for (Int_t m=0; m<fgMaxInput; m++)
1406   {
1407     if (measured->GetBinContent(m+1) < 1)
1408       continue;
1409       
1410     for (Int_t loop=0; loop<4; loop++)
1411     {
1412       //Printf("%d %d", i, loop);
1413       Float_t factor = 1;
1414       switch (loop)
1415       {
1416         case 0: factor = 0.5; break;
1417         case 1: factor = -0.5; break;
1418         case 2: factor = 1; break;
1419         case 3: factor = -1; break;
1420         default: return -1;
1421       }
1422       
1423       TH1* measuredClone = (TH1*) measured->Clone("measuredClone");
1424   
1425       measuredClone->SetBinContent(m+1, measured->GetBinContent(m+1) + factor * TMath::Sqrt(measured->GetBinContent(m+1)));
1426       //new TCanvas; measuredClone->Draw("COLZ");
1427     
1428       newResult[loop]->Reset();
1429       Unfold(correlation, efficiency, measuredClone, measuredClone, newResult[loop]);
1430       // WARNING if we leave here, then nothing is calculated
1431         //return -1;
1432       
1433       delete measuredClone;
1434     }
1435     
1436     for (Int_t t=0; t<fgMaxParams; t++)
1437     {
1438       // one value
1439       //matrix(i, j) = (result->GetBinContent(j+1) - newResult->GetBinContent(j+1)) / TMath::Sqrt(changedHist->GetBinContent(1, i+1));
1440       
1441       // four values from the derivate (procedure taken from ROOT -- suggestion by Ruben)
1442       // = 1/2D * [ 8 (f(D/2) - f(-D/2)) - (f(D)-f(-D)) ]/3
1443       
1444       /*
1445       // check formula
1446       measured->SetBinContent(1, m+1, 100);
1447       newResult[0]->SetBinContent(t+1, 5 * 0.5 + 10);
1448       newResult[1]->SetBinContent(t+1, 5 * -0.5 + 10);
1449       newResult[2]->SetBinContent(t+1, 5 * 1 + 10);
1450       newResult[3]->SetBinContent(t+1, 5 * -1 + 10);
1451       */
1452       
1453       matrix(m, t) = 0.5 / TMath::Sqrt(measured->GetBinContent(m+1)) * 
1454         ( 8. * (newResult[0]->GetBinContent(t+1) - newResult[1]->GetBinContent(t+1)) - 
1455               (newResult[2]->GetBinContent(t+1) - newResult[3]->GetBinContent(t+1))
1456         ) / 3;
1457     }
1458   }
1459   
1460   for (Int_t loop=0; loop<4; loop++)
1461     delete newResult[loop];
1462   
1463   //matrix.Print();
1464   //TH2* matrixHist = new TH2D(matrix); new TCanvas; matrixHist->Draw("BOX");
1465   
1466   // ...calculate folded guess  
1467   TH1* convoluted = (TH1*) measured->Clone("convoluted");
1468   convoluted->Reset();
1469   for (Int_t m=0; m<fgMaxInput; m++)
1470   {
1471     Float_t value = 0;
1472     for (Int_t t = 0; t<fgMaxParams; t++)
1473     {
1474       Float_t tmp = correlation->GetBinContent(t+1, m+1) * result->GetBinContent(t+1);
1475       if (efficiency)
1476         tmp *= efficiency->GetBinContent(t+1);
1477       value += tmp;
1478     }
1479     convoluted->SetBinContent(m+1, value);
1480   }
1481   
1482   // rescale
1483   convoluted->Scale(measured->Integral() / convoluted->Integral());
1484   
1485   //new TCanvas; convoluted->Draw(); measured->Draw("SAME"); measured->SetLineColor(2);
1486   //return;
1487   
1488   // difference
1489   convoluted->Add(measured, -1);
1490   
1491   // ...and the bias
1492   for (Int_t t = 0; t<fgMaxParams; t++)
1493   {
1494     Double_t error = 0;
1495     for (Int_t m=0; m<fgMaxInput; m++)
1496       error += matrix(m, t) * convoluted->GetBinContent(m+1); 
1497     result->SetBinError(t+1, error);
1498   }
1499   
1500   //new TCanvas; result->Draw(); gPad->SetLogy();
1501   
1502   return 0;
1503 }